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ABSTRACT 



Aims. The Earth's atmosphere affects ground-based astronomical observations. Scattering, absoiption, and radiation processes dete- 
riorate the signal-to-noise ratio of the data received. For scheduling astronomical observations it is, therefore, important to accurately 
estimate the wavelength-dependent effect of the Earth's atmosphere on the observed flux. 

Methods. In order to increase the accuracy of the exposure time calculator of the European Southern Observatory's (ESO) Very Large 
Telescope (VLT) at Cerro Paranal, an atmospheric model was developed as part of the Austrian ESO In-Kind contribution. It includes 
all relevant components, such as scattered moonlight, scattered starlight, zodiacal light, atmospheric thermal radiation and absorption, 
and non-thermal airglow emission. This paper focuses on atmospheric scattering processes that mostly affect the blue (< 0.55 /jm) 
wavelength regime, and airglow emission lines and continuum that dominate the red (> 0.55 fim) wavelength regime. While the 
former is mainly investigated by means of radiative transfer models, the intensity and variability of the latter is studied with a sample 
of 1 186 VLT FORS 1 spectra. 

Results. For a set of parameters such as the object altitude angle, Moon-object angular distance, ecliptic latitude, bimonthly period, 
and solar radio flux, our model predicts atmospheric radiation and transmission at a requested resolution. A comparison of our model 
with the FORS 1 spectra and photometric data for the night-sky brightness from the literature, suggest a model accuracy of about 20%. 
This is a significant improvement with respect to existing predictive atmospheric models for astronomical exposure time calculators. 
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1. Introduction 

Ground-based astronomical observations are affected by the 
Earth's atmosphere. Light from astronomical objects is scattered 
and absorbed by air molecules and aerosols. This extinction ef- 
fect can cause a significant loss of flux, depending on the wave- 
length and weather conditions. The signal of the targeted object 
is further deteriorated by background radiation, which is caused 
by light from other astronomical radiation sources scattered into 
the line of sight and emission originating from the atmosphere 
itself. Since these contributions can vary significantly with time, 
the achievable signal-to-noise ratio for an astronomical observa- 
tion strongly depends on the state of the Earth's atmosphere and 
the Sun-Earth-Moon system. Therefore, for efficient time man- 
agement of any modern observatory, it is critical to provide a 
reliable model of the Earth's atmosphere for estimating the expo- 
sure time required to achieve the goals of scientific programmes. 
Data calibration and reduction also benefit from a good knowl- 
edge of atmospheric effects (see e.g. Davies l2~007l) . 

For this reason, various investigations were performed to 
characterise the atmospheric conditions at telescope sites (see 
Leinert et al. 1998 for a comprehensive overview). Photometric 
measurements of the night-sky brightness and its variability 
were done at e.g. Mauna Kea (Krisciunas l 19971) . La Palma (Benn 
& Ellison [1998) , Cerro Tololo (Walker [T9871 Krisciunas et al. 



* Based on observations made with ESO telescopes at Paranal 
Observatory 



120071 La S illa (Mattila et al. [19961 1, a nd Cer ro Paranal (Patat 
120031 120081 . For Cerro Paranal, Patat (120081 also carried out 
a detailed spectroscopic analysis, and found that the night sky 
showed strong variations of more than one magnitude. Also, the 
sky brightness depends on the solar activity cycle (Walker 1988 1 
Patat l2008l l. It is related to the variations of the upper atmosphere 
airglow line and continuum emission, which dominate the near- 
UV, optical, and near-IR sky emission under dark-sky conditions 
(Chamberlain [T9r5T1 R oach & Gordon [T973l Leinert et al. [T998l 
Khomich et al. 120081) . When the Moon is above the horizon, 
scattered moonlight dominates the blue wavelengths (Krisciunas 
& Schaefer |1991| l. A much weaker, but always present, compo- 
nent is scattered starlight. The distribution of integrated starlight 
(Mattila [T980l Toller[T9"8"Tl Toller et al. [T987l Leinert et al.[1998 1 
Melchior et al. 2007 ) and how it is scattered in the Earth's atmo- 
sphere has been studied (Wolstencroft & van Breda 1967 Staude 
119751 ; Bernstein et al. |2002). A significant component at opti- 
cal wavelengths is the so-called zodiacal light, solar radiation 
scattered by interplanetary dust grains mainly distributed in the 
ecliptic plane (Levasseur-Regourd & Dumont 1980 Mattila et 
al. [1996; Leinert et al. 1998). For a realistic description of the 
zodiacal light intensity distribution for ground-based observa- 
tions, scattering calculations are also required (Wolstencroft & 
van Breda [19671 Staude [19751 Bernstein et al. 120021 . Finally, 
the wavelength-dependent extinction of radiation from astro- 
nomical objects by Rayleigh scattering off of air molecules, Mie 
scattering off of aerosols, and absorption by tropospheric and 
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Table 1. Sky model parameters for optical wavelength range 



Parameter 


Description 


Unit 


Range 


Default* 


Demo run'' 


Section 


90° - zo 


altitude of target above horizon 


deg 


[0,90] 


90. 


85.1 




a 


separation of Sun and Moon as seen from Earth 


deg 


[0,180] 


0. 


77.9 


im 


o 
H 


separation of Moon and target 




[0,180] 


180. 


51.3 


i3~n 


90 Zmoon 


altitude of Moon above horizon 


deg 


[-90,90] 


-90. 


41.3 


ED 


^raoon 


relative distance to Moon (mean =1) 




[0.945,1.055] 


1. 


L d 


ED 


A-Aq 


heliocentric ecliptic longitude of target 


deg 


[-180,180] 


135. 


-124.5 


1331 


P 


ecliptic latitude of target 


deg 


[-90,90] 


90. 


-31.6 


[33] 


S 10.7 cm 


monthly- averaged solar radio flux at 10.7 cm 


sfu c 


> 


130. 


205.5 


14721-1441 


P 

1 season 


bimonthly period (1: Dec/Jan, 6: Oct/Nov; 0: entire year) 




[0,6] 





4 


I2~f1l473ll44l 


p 

1 lime 


period of the night (x/3 of night, x = 1,2,3; 0: entire night) 




[0,3] 





3 


14311441 


vac/air 


vacuum or air wavelengths 




vac/air 


vac 


air 


im 1431 1441 



Notes. 

(fl) We neglect temperature and emissivity of telescope and instrument, because these parameters are irrelevant for the optical spectral range. 

{b) Used for Tabled] 

(c) Used for Figs. 0|6] [13] andQjJ 

{d} Fixed to default value because of its minor importance (see also Sect. l3.lt . 
w 1 sfu = 0.01 MJy. 




1 2 3 

Wavelength [/um] 



Fig. 1. Components of the sky model in logarithmic radiance 
units for wavelengths between 0.3 and 4.2yum. This example, 
with the Moon above the horizon, shows scattered moonlight, 
scattered starlight, zodiacal light, thermal emission by the tele- 
scope and instrument, molecular emission of the lower atmo- 
sphere, airglow emission lines of the upper atmosphere, and air- 
glow/residual continuum. The scattered light and airglow com- 
ponents were only computed up to the K band because of their 
negligible importance at longer wavelengths. For the model pa- 
rameters used, see Table [TJ 



stratospheric molecules was studied and characterised at differ- 
ent telescope sites, such as La Silla (Tug 1977, 1980, Rufener 
[19861 Sterken & Ma nfroid[T992l Burki et al. 2005), Cerro Tololo 
(Stone & Ba ldwin [19831 Baldwin & Stone [19841 Gutierrez- 
Moreno et al. rT982irT986l >, and Cerro Paranal (Patat et al. l20TT1 i. 

Due to the complexity and variability of the night-sky radi- 
ation, a good atmospheric radiation model is crucial for a re- 
liable astronomical exposure time calculator (ETC). Currently, 
the European Southern Observatory (ESO) uses a sky back- 
ground model for its ETC (Ballester et al. 2000 1, which includes 
the photometric night-sky brightness measurements of Walker 



([19871) for the optical and Cuby et al. ([2000) for the near-IR. The 
U to I magnitudes of Walker vary as a function of Moon phase. 
For optical spectrographs, a night-sky spectrum is calculated by 
scaling an observed, instrument-dependent template spectrum 
to the Walker filter fluxes. The spectroscopic near-IR/mid-IR 
model is based on the line atlas of Hanuschik (2003) and the 
OH airglow calculations of Rousselot et al. (2000), thermal tele- 
scope emission, an instrument-related constant continuum, and 
atmospheric thermal line and continuum emission. The latter is 
provided as a set of template spectra computed with the radiative 
transfer code Reference Forward Modefjfor different airmasses 
and water vapour column densities. Since the current ETC is 
limited in reproducing the variable intensity of the night sky, 
we have developed an advanced model for atmospheric radiation 
and transmission, which includes scattered moonlight, scattered 
starlight, zodiacal light, thermal emission from the telescope, 
molecular emission and absorption in the lower atmosphere, and 
airglow line and continuum emission, including their variability 
with time. This "sky model" has been derived for Cerro Paranal 
in Chile (2635 m, 24° 38' S, 70° 24' W). It is also expected 
to work for the nearby Cerro Armazones (3064m, 24° 36' S, 
70° 12' W), the future site of the European Extremely Large 
Telescope (E-ELT), without major adjustments. An example of 
a spectrum computed by our sky model is given in Fig. [TJ The 
crucial input parameters are listed in Table [TJ The model will 
be incorporated into the ESO ETC0 and will be made publicly 
available to the community via the ESO web site. 

This paper focuses on the discussion of model components 
relevant for the wavelength range from 0.3 to 0.92 /im, which is 
designated as "optical" in the following. We will discuss scat- 
tering and absorption in the atmosphere (Sect. [2}, contributions 
from extraterrestrial radiation sources (Sect. 0], and the inten- 
sity and variability of airglow emission lines and continuum 
(Sect. @J. The quality of our optical sky model will be evalu- 
ated in Sect. [5] The near-IR and mid-IR regimes will be treated 
in a subsequent paper. 



1 http : // www . atm . ox . ac . uk/RFM/ 

2 http : //www . eso . org/observing/etc/ 
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If intensity units are not explicitly given in this pa- 
per, photons s~ m _2 jUm _1 arcsec~ 2 are taken as the standard. 
Magnitudes are always given in the Vega system. 



2. Atmospheric extinction 

Light from astronomical objects is scattered and absorbed in the 
Earth's atmosphere. For point sources, both effects result in a 
loss of radiation, which is usually described by a wavelength- 
dependent extinction curve or an atmospheric transmission 
curve. The components of such a curve for Cerro Paranal are 
discussed in Sect. 12.11 For extended sources, light is not only 
scattered out of the line of sight, it is also scattered into it. This 
effect causes an effective extinction curve to differ from that of 
a point source and depends on the spatial distribution of the ex- 
tended emission. To quantify the change of the extinction curve, 
we performed three-dimensional (3D) scattering calculations, 
which are described in Sect. 12.21 

2.1. The transmission curve 

The atmospheric transmission depends on scattering and absorp- 
tion. Light can be scattered in dry atmosphere by air molecules 
such as N2 and O2 or by aerosols like silicate dust, sea salt, 
soot, or droplets of sulphuric acid. The former effect is known 
as Rayleigh scattering, and is characterised by a strong wave- 
length dependence proportional to A ~ 4 and a relatively isotropic 
scattering phase function (a factor of 2 variation for unpolarised 
radiation). The latter can be described by Mie scattering if spher- 
ical particles are assumed. Aerosol scattering is characterised by 
a relatively weak wavelength dependence (A~ l to /L 2 ) and pro- 
nounced forward scattering, for which the maximum intensity 
can easily be two orders of magnitude higher than at large scat- 
tering angles. Atmospheric absorption in the optical is mainly 
caused by bands from three molecules: molecular oxygen (O2), 
water vapour (H2O), and ozone (O3). The absorption of O2 is a 
function of atmospheric density. While water absorption is most 
efficient close to the ground, where the absolute humidity is the 
largest, ozone mainly absorbs at stratospheric altitudes of about 
20 km. The combination of scattering and absorption results in 
an extinction curve, which is often given in mag airmass -1 , or a 
transmission curve providing values from (totally opaque) to 
1 (fully transparent). The transmission t(A) can be linked to the 
zenithal optical depth tq(A) and zenithal extinction coefficient 
k(A) by 

t(A) = e- T °»* = l -°- 4 *« x . (1) 

The airmass X can be calculated by the formula of Rozenberg 
(fl966l l: 

X = (cos(z) + 0.025 e~ 1 1 cosfc) )~ ' , (2) 

where z is the zenith distance and X converges to 40 at the hori- 
zon. 

For Cerro Paranal, Fig. |2] shows the annual mean transmis- 
sion curve at zenith and its components. The extinction at blue 
wavelengths is dominated by Rayleigh scattering. This compo- 
nent is very stable and can be well reproduced by the parametri- 
sation 



tr(A) = 



1013.25 



(0.00864 + 6.5 x 10~ 6 //) 



x A 



-(3.916 + 0.074 A + 0.050 / A) 



(3) 



a 
to 

H 




0.6 0.8 1 

Wavelength [/itm] 

Fig. 2. Annual mean zenith transmission curve for Cerro Paranal 
(black solid line). The Earth's atmosphere extinguishes the flux 
from point sources by Rayleigh scattering by air molecules (red 
dash-dotted line), Mie scattering by aerosols (green dashed line), 
and molecular absorption (blue solid line). For the plotted wave- 
length range, the latter is caused by molecular oxygen (A band 
at 0.762 fim, B band at 0.688 /mi, and y band at 0.628 fim), water 
vapour (prominent bands at 0.72, 0.82, and 0.94 /urn), and ozone 
(Huggins bands in the near-UV and broad Chappuis bands at 
about 0.6 fim). 




Feb/Mar average 
Aug/Sep average 
Annual 1(7 dev. 

_] I I I L_ 



with wavelength A in //m (see Liou [20021 . For the pressure p 
and the height H, we take 744 hPa and 2.64km respectively. 



0.8 0.9 1 

Wavelength [/iim] 

Fig. 3. Variation of molecular absorption for Cerro Paranal. The 
extreme bimonthly mean transmission curves and 1 <x deviations 
of the annual mean curve (red outer curves) are shown. The high- 
est mean transmission (lowest water vapour content) is found 
for August/September (green curve), and the lowest arises in 
February /March (black curve). In contrast to H2O, the variations 
of the O2 bands are very small. For the identity of the bands, see 
Fig.E] 



The pressure corresponds to the annual mean for Cerro Paranal 
(743.5 ±1.5 hPa), as derived from the meteorological station of 
the VLT Astronomical Site Monitor. 

At red wavelengths, aerosol scattering becomes as impor- 
tant as Rayleigh scattering. However, the total amount of ex- 
tinction by scattering is small in this wavelength regime. For 
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Cerro Paranal, Patat et al. (1201 1) provide an approximation for 
the aerosol extinction derived from 600 VLT FORS 1 spectra 
observed over six months. The aerosol extinction coefficient is 
parametrised by 

kM«k A a , (4) 

where k = 0.013 + 0.002 mag airmass" 1 and a = -1.38 + 0.06, 
with the wavelength A in jjm. Due to an increased discrepancy 
between the fit and the observed data in the near-UV (see Patat 
et al. 1201 11 1, we use Eq. |4] only for wavelengths longer than 
0.4 fim. For shorter wavelengths, we use a constant value of 
^aer = 0.050 mag airmass -1 , which corresponds to the fit value at 
0.4 fim. The density, distribution, and composition of aerosols is 
much more variable than what is observed for the air molecules, 
which determine the Rayleigh scattering. Patat et al. (201 find 
that k- der varies by about 20% at 0.4 fim. However, these varia- 
tions are of minor importance for the total transmission of the 
Earth's atmosphere, since the aerosol extinction coefficients are 
small. 

Figure [2] exhibits several prominent absorption bands (see 
also Patat et al. 1201 II) . At wavelengths below 0.34 fim, there is 
a conspicuous fall-off of the transmission curve caused by the 
Huggins bands of ozone. The stratospheric ozone layer is also 
responsable for the Chappuis absorption bands between 0.5 and 
0.7 fim. The relatively narrow, but strong, bands at 0.688 fim and 
0.762 fim can be identified as the Fraunhofer B and A bands of 
molecular oxygen. Finally, the complex bands at 0.72, 0.82, and 
0.94 fim are produced by water vapour. 

The molecular absorption bands have been calculated using 
the Line By Line Radiative Transfer Model (LBLRTM), an at- 
mospheric radiative transfer code provided by the Atmospheric 
and Environmental Research Inc. (see Clough et al. 2005). 
This widely used code in the atmospheric sciences computes 
transmission and radiance spectra based on the molecular line 
database HITRAN (see Rothman et al. 120091) and atmospheric 
vertical profiles of pressure, temperature, and abundances of rel- 
evant molecules. 

For Cerro Paranal, we use merged atmospheric profiles from 
three data sources to reproduce the climate and weather con- 
ditions in an optimal way. First, the equatorial daytime stan- 
dard profile from the MIPAS instrument of the ENVIS AT sat el- 
lite (prepared by J. Remedios 2001; see Seifahrt et al. 12010) is 
taken. It provides abundances for 30 molecular species up to 
an altitude of 120 km. Following Patat et al. (1201 11 ). the ozone 
profile is corrected by a factor of 1.08 to achieve a column 
density of 258 Dobson units, which represents the mean value 
for Cerro Paranal. Second, we use profiles from the Global 
Data Assimilation SysterrQ (GDAS), maintained by the Air 
Resources Laboratory of the National Oceanic and Atmospheric 
Administration (cf. Seifahrt et al. 120101 The GDAS profiles for 
pressure, temperature, and relative humidity are provided on a 
3 h basis for a 1 ° x 1 global grid. These models are adapted to 
data from weather stations all over the world and are suitable for 
weather-dependent temperature and water vapour profiles up to 
altitudes of 26 km. Third, data from the meteorological station 
at Cerro Paranal are used to scale the pressure, temperature, and 
water vapour profiles at the altitude of the mountain. For higher 
altitudes, the scaling factor is reduced and approaches 1 at 5 km. 

For our sky model, we analysed the resulting data set and 
constructed mean profiles with their 1 <x deviations for different 
periods. We divide the year into six two-month periods, start- 
ing with December/January (cf. Table Q] and Sect. [4](. The two 
most extreme mean spectra are shown in Fig. [3] The seasonal 

3 http : //ready . ar 1 . noaa . gov/gdas 1 . php 



Observer's 




Top of 
atmosphere 



Fig. 4. Geometry of the scattering in the Earth's atmosphere (cf. 
Wolstencroft & van Breda |T967j ). Point N is at the top of atmo- 
sphere and not in the same plane as the other points. The azimuth 
of N as seen from S is A (not shown). S and T are at an azimuth 
Aq (also not shown) for an observer at O. 



variability of the H2O bands is clearly visible. For the driest pe- 
riod (August/September), the mean absorption is only half the 
amount of the most humid period (February /March). The total 
intra-annual variability indicates line depth variations of an order 
of magnitude, i.e. large statistical uncertainties. For this reason, 
the significant seasonal dependence is included in the sky model 
and the less pronounced average, nocturnal variations have been 
neglected. Apart from the two-month period, only the airmass 
(see Eq. [2]) is used as input parameter for the computation of 
transmission curves. Since radiative transfer calculations with 
LBLRTM are time-consuming, the sky model is run with a pre- 
calculated library of transmission spectra, consisting of spectra 
for the different bimonthly periods and a regular grid of five air- 
masses between 1 and 3. This is sufficient for a reliable interpo- 
lation of the airmass-dependent change of spectral features. 

A more detailed discussion on atmospheric profiles, radiative 
transfer codes, and the properties of transmission and radiance 
spectra, especially in the near-IR and mid-IR, will be given in a 
subsequent paper. 

2.2. Atmospheric scattering 

The estimate of scattered light from extended sources, such as 
integrated starlight, zodiacal light, or airglow in the Earth's at- 
mosphere, requires radiative transfer calculations. Since the op- 
tical depth for scattering is relatively small in most of the opti- 
cal wavelength range (see Fig. single-scattering calculations 
provide a sufficient approximation. In this case, the computa- 
tions can be performed in 3D with a relatively compact code 
(see Wolstencroft & van Breda 1 19671 Staude |TW5l Bernstein et 
al. l2002T ). 
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To obtain the integrated scattered light towards the azimuth 
Aq and zenith distance zo, we consider scattering path elements 
S of density n(cr), with cr being the radius vector from the centre 
of Earth C to S , from the top of atmosphere T to the observer O 
at height Ho above the surface (see Fig.|4]i. The distance of O to 
C is ctq = Hq + R, where R is the radius of Earth (= 6371 km for 
the mean radius). For each path element S at distance S2 from O, 
the contributions of radiation from all directions (A, z) to the in- 
tensity at S are considered. The integration over the solid angle 
depends on the spatial intensity distribution Iq(A,z) of the ex- 
tended radiation source and the path of the light from the entry 
point in the atmosphere N to S of length s\ . The latter includes 
the scattering of light out of the path (and possible absorption), 
which depends on the effective column density of the scatter- 
ing/absorbing particles 



transfer calculations in a plane-parallel atmosphere (see Dave 
[19641 Staude [l975l We multiply 7 scat by the factor 



Bi(z, a) 



-I 



site cr) 



n(cr') dcr'. 



(5) 



It also depends on the wavelength-dependent extinction cross 
section C ext (A), which can be derived from the optical depth at 
zenith To (see Eq.Q]) by means of Eq. [5] and 



C ext (A) - 



tq(A) 
Bo(0,cr y 



(6) 



After the scattering at S , the intensity is further reduced by the 
effective column density ,62(20, cr) along the path Si. Thus, the 
calculations can be summarised by 

/scat(A ,Zo) = ^ I I I n(cr)P(0) 

x/ (A,z)e" c °' ( ' ,)[iilfcfr)+S2(a ' <r)] dA sinzdzds. (7) 

C SC at is the wavelength-dependent scattering cross section, which 
will deviate from C ex t if absorption occurs. The maximum zenith 
distance z max is higher than 90° and depends on the height of S 
above the ground. The scattering phase function P depends on 
the scattering angle 9, the angle between the paths S\ and 

Rayleigh scattering (see Sect. 12. U is characterised by the 
phase function 

P(0) = !(l+cos 2 (0)). (8) 



Similar to Bernstein et al. (120021 1. we neglect the effect of polar- 
isation on the scattering phase function. Even for zodiacal light, 
where some polarisation is expected (e.g. Staude IT9751 I. the po- 
larisation does not appear to significantly affect the integrated 
scattered light. Results of Wolstencroft & van Breda (1967 1 sug- 
gest deviations of only a few per cent. For the vertical distribu- 
tion of the scattering molecules, we use the standard barometric 
formula 

n(h) = noe- h/b( >. (9) 

Here, h — cr' — R, the sea level density «o = 2.67 x 10 19 cm" 3 , 
and the scale height ho - 7.99 km above the Earth's surface (cf. 
Bernstein et al. 2002). For the troposphere and the lower strato- 
sphere, where most of the scattering occurs, this is a good ap- 
proximation. Cerro Paranal is at an altitude of Hq = 2.64 km. 
For the thickness of the atmosphere, we take H max = 200 km. 
With these values, C ext is 1.75 x 10~ 26 cm 2 for t = 0.27, which 
corresponds to a wavelength of about 0.4 jjm. For Rayleigh scat- 
tering, C SC at = Cext, i-e. no absorption is involved. In the near-UV, 
where the zenithal optical depth is greater than 0.3, the single 
scattering approximation becomes questionable. For this reason, 
we apply a multiple scattering correction derived from radiative 



'MS 



1 +2.2t . 



(10) 



The uncertainty is in the order of 5%. The factor Fms does not in- 
clude reflection from the ground. Mattila (2003 1 reports a ground 
reflectance of about 8% in the region near the Las Campanas 
Observatory in autumn. Since this telescope site is also located 
in the Atacama desert, we assume this value with an uncertainty 
of a few per cent. Using the tables of Ashburn (11954) , an 8% 
reflectance translates into an additional 7 scal correction factor of 
about 1.07. 

The Mie scattering of aerosols is characterised by a phase 
function with a strong peak in forward direction. For this study, 
we take the measured P(6) of Green et al. (|1971| l, which cov- 
ers the peak up to an angle of 20°. The phase function for larger 
scattering angles is extrapolated by a linear fit in the log P - log 9 
diagram. This simple approach neglects the increase of the phase 
function for scattering angles close to 180° (see e.g. Liou|2002). 
However, since P already decreases by a factor of 30 from 0° to 
20°, the details of the phase function at larger angles are not cru- 
cial for the scattering of extended emission. The total scattered 
intensity is completely dominated by the contribution from an- 
gles close to the forward direction. For the height distribution of 
aerosols, we also take Eq.|9] but assume Hq = 1 . 1 1 x 10 4 crrT 3 and 
ho = 1.2 km. This is the tropospherical distribution of Elterman 
( 1 19661 ). Following Staude ( 119751 ), we neglect the stratospheric 
aerosol component, which contains only about 1% of the par- 
ticles. Using Eq. |6] we obtain C ext = 3.39 x 10~ 10 cm 2 for 
to = 0.05, which characterises wavelengths of about 0.4 /mi. 
Dust, in particular soot, also absorbs radiation. For this reason, 
C S cat is lower than C ex t- The ratio of these is called the single 
scattering albedo at. Following the recommendation of Mattila 
(2003), ci) = 0.94 is used for the model, i.e. strong absorbers like 
soot do not significantly contribute to the aerosol population. We 
neglect any corrections for multiple scattering and ground re- 
flection for Mie scattering, since the optical depths are low at all 
relevant wavelengths (see Fig. |2]i and forward scattering domi- 
nates. 

To test our model, we computed scattering intensities for a 
uniform radiation source of unit brightness. For the surface level 
and To = 0.052, we obtain scattering fractions of 0.024 for the 
zenith and 0.360 for the horizon (zo = 90°) with the multiple 
scattering corrections neglected. For Rayleigh scattering, these 
results can be compared to those of Staude (119751 I. Our values 
are only about 5% higher. This is a good agreement compared 
to differences between Staude (119751 1 and older calculations of 
Wolstencroft & van Breda (119671 1 and Ashburn (11954) . Figure |5] 
shows the resulting scattering fractions for Rayleigh and Mie 
scattering for different optical depths, which have been con- 
verted into wavelengths based on the Cerro Paranal mean trans- 
mission curve discussed in Sect. 12. II For zenith distances that are 
realistic for astronomical observations, wavelengths longer than 
about 0.45 fim for Rayleigh scattering, and all wavelengths for 
aerosol scattering, the scattering fractions stay below 10% for 
a source of uniform brightness. Scattering by air molecules and 
the ground surface can be important for lines of sight close to the 
horizon and near-UV wavelengths. Figure ^indicates scattering 
fractions of about 40% under these conditions. Moreover, there 
appears to be a downturn of the scattering fraction for very high 
zenith distances (see also Staude 1975), where the line-of-sight 
direction becomes nearly opaque. In this case, the light that is 
not extinguished has to be scattered near the observer and the 
incident photons have to originate close to the zenith. 
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Fig. 5. Intensity of scattered radiation for a uniform source of 
unit brightness using the Paranal extinction curve (see Sect. 12. U . 
Contributions of scattered light into the line of sight are shown 
for Rayleigh (solid lines) and Mie scattering (dashed lines) and 
zenith distances from 0° to 80° in 10° steps plus an extreme value 
of 85°. Except for Rayleigh scattering at very short wavelengths, 
the scattering intensities increase with zenith distance. 



The results of scattering calculations for more realistic 
source intensity distributions will be discussed in Sects. [3T2l 13 . 31 
andPO 



3. Moon, stars, and interplanetary dust 

Bright astronomical objects affect night-sky observations by the 
scattering of their radiation into the line of sight. The relevant 
initial radiation sources are the Sun and the summed-up light 
of all other stars. The solar radiation must first be scattered by 
interplanetary dust grains and the Moon surface to contribute 
to the night-sky brightness. In the following, we discuss these 
scattering-related sky model components, i.e. scattered moon- 
light (Sect. lXTT ). scattered starlight (Sect. [372), and zodiacal light 
(Sect. [331 ). Another source of scattered light is man-made light 
pollution, which can be neglected for Cerro Paranal. Patat (2003) 
did not detect spectroscopic signatures of such contamination. 

3.1. Scattered moonlight 

For observing dates close to Full Moon, scattered moonlight is 
by far the brightest component of the optical night-sky back- 
ground. In particular, blue wavelengths are affected, since the 
Moon spectrum resembles a solar spectrum and the efficency of 
Rayleigh/Mie scattering increases towards shorter wavelengths 
(see Fig.UJi. 

The Moon can be considered a point source. For this rea- 
son, it is not necessary to perform the scattering calculations 
for extended sources, as discussed in Sect. I2.2I Instead, a semi- 
analytical model for the photometric V band by Krisciunas & 
Schaefer (1991) is used for the sky model, which has been ex- 
tended into a spectroscopic version. 

Based on Krisciunas & Schaefer ( 119911 ), we compute the 
Moon-related sky surface brightness as the sum of the contri- 
butions from Rayleigh and Mie scattering 
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Fig. 6. Effects of atmospheric scattering on radiation originating 
outside of the atmosphere. The displayed radiation sources are 
Moon (blue solid lines), stars (green dashed lines), and interplan- 
etary dust (red dash-dotted lines). All spectra are normalised to 
unity at 0.5 pm. In the upper panel, a solar spectrum is assumed 
as top-of-atmosphere spectrum for moonlight (see Sect. I3.lt . 
A slightly reddened solar spectrum is taken for solar radiation 
scattered at interplanetary dust grains (see Sect. 13.3) . Integrated 
starlight tends to be redder than sunlight at long wavelengths (see 
Sect. 13. 2t . The lower panel shows the resulting low-resolution 
spectra, after light has been scattered in the Earth's atmosphere, 
for observing conditions as listed in Table [1] While the lunar and 
stellar contributions represent scattered radiation only, the dust- 
related zodiacal light consists of scattered and direct light. 



»(T)(i-f*; M M)). 



(12) 

(13) 
(14) 



Bmoon(^) - #moon,R(T) + B„ 



(ID 



where 

Bmoon.R/M^) = /r/m(p) 

The empirical scattering functions 

/ R (p)= 10 5 - 7() (l.06 + cos 2 (p)) 

for Rayleigh scattering and 

/ M (p) = 10 7 - 15 -^ 

for Mie scattering depend on the angular separation of Moon 
and object p, which is restricted to angles greater than 10°. The 
functions deviate from the ones of Krisciunas & Schaefer (119911 1 
by factors of 2.2 and 10 respectively. These corrections are nec- 
essary because of the separation of the two scattering processes 
in Eqs. Q~T|and[T2] the model extension in wavelength, and an 
optimisation for the observing conditions at Cerro Paranal. The 
two correction factors were derived separately from a compari- 
son with optical VLT data (see Sect. 14. 2t by using the increasing 
importance of Rayleigh scattering with respect to Mie scatter- 
ing for shorter wavelengths and larger scattering angles p (see 
Sect. I2.21 >. The factor of 10 for Mie scattering is highly uncertain 
due to a lack of clear constraints from the available data set. The 
Moon illuminance is proportional to 

/* oc 10 - 0.4 (0.026 10| + 4.0 X10-V 4 ) x 



-*moon 



i) - 



(15) 
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where the Moon distance d moon can vary up to 5.5% relative to 
the mean value. The lunar phase angle a = 180° - is the sep- 
aration angle of Moon and Sun as seen from Earth. It can be 
estimated from the fractional lunar illumination (FLI) by 

a « arccos(l -2FLI), (16) 

since the influence of the lunar ecliptic latitude on a is negligi- 
ble. B moon also depends on the atmospheric transmission for the 
airmass of the Moon t Xmmn , for which the Patat et al. (1201 11 1 ex- 
tinction curve is used (see Fig. |2). Lastly, it depends on a term 
that describes the amount of scattered moonlight in the viewing 
direction, which can be roughly approximated by 1 - t*° M for 
either Rayleigh or Mie scattering and the target airmass Xq. The 
airmasses of the Rrisciunas & Schaefer (119911 ) model are com- 
puted using 

X = (l - 0.96 sin 2 (z))~°' 5 , (17) 

which is especially suited for scattered light, because of a low 
limiting airmass at the horizon (cf. Eq.[2). 

The Krisciunas & Schaefer model was developed only for 
the V band. It can be extended to other wavelength ranges by 
using the Cerro Paranal extinction curve (see Sect. 12. 11 and its 
Rayleigh and Mie components for the transmission curves t, ?r, 
and ?m in Eq. Q~2] and by assuming that the Moon spectrum re- 
sembles the spectrum of the Sun (see Fig.[6](. The latter is mod- 
elled by the solar spectrum of Colina et al. ( 1996) scaled to the 
V-band brightness of the moonlight scattering model. For most 
situations, this approach results in systematic uncertainties of no 
more than 10% of the total night-sky flux compared to observa- 
tions. The statistical uncertainties are in the order of 10 to 20% 
(see Sect. [5}. The errors become larger for Full Moon and/or 
target positions close to the Moon (p % 30°). The latter espe- 
cially concerns red wavelengths, where the contribution of Mie 
forward scattering to the total intensity of scattered light is par- 
ticularly high. Consequently, the strongly varying aerosol prop- 
erties have to be known (which is a challenging task) to allow 
for a good agreement of model and observations. However, op- 
tical astronomical observations close to the Moon are unlikely. 
Therefore, this is not a critical issue for the sky model. 

3.2. Scattered starlight 

Like moonlight, starlight is also scattered in the Earth's atmo- 
sphere. However, stars are distributed over the entire sky with a 
distribution maximum towards the centre of the Milky Way. This 
distribution requires the use of the scattering model for extended 
sources described in Sect. 12.21 Since scattered starlight is only a 
minor component compared to scattered moonlight and zodiacal 
light (see Fig. |T}, it is sufficient for an ETC application to com- 
pute a mean spectrum. This simplification allows us to avoid the 
introduction of sky model parameters such as the galactic coor- 
dinates, which have a very low impact on the total model flux. 

For the integrated starlight (ISL), we use Pioneer 10 data at 
0.44 jt/m (Toller [19811 Toller et al. [T9871 Leinert et al. [1998). 
These data are almost unaffected by zodiacal light. The small 
"hole" in the data set towards the Sun has been filled by inter- 
polation. Since the Pioneer 10 maps do not include stars brighter 
than 6.5 mag in V, a global correction given by Melchior et al. 
(2007) is applied, which increases the total flux by about 16%. 
Scattering calculations are only performed for the distribution of 
starlight in the B band, i.e. the dependence of this distribution 
on wavelength is neglected. This approach is supported by cal- 
culations of Bernstein et al. (2002), indicating that the shapes 



of ISL spectra are very stable, except for regions in the dusti- 
est parts of the Milky Way plane. The wavelength dependence 
of the scattered starlight is considered by multiplying the rep- 
resentative ISL mean spectrum of Mattila (1980) by the result- 
ing wavelength-dependent amount of scattered light (for illustra- 
tion see Figs. [5] and [6]). Since the mean ISL spectrum of Mattila 
only covers wavelengths up to 1 fan, we extrapolated the spec- 
tral range by fitting an average of typical spectral energy distri- 
butions (SED) for early- and late-type galaxies produced by the 
SED-fitting code CIGALE (see Noll et al. l2009b . Since the SED 
slopes are very similar in the near-IR, the choice of the spectral 
type is not crucial. The solar and the mean ISL spectrum are sim- 
ilar at blue wavelengths, but the ISL spectrum has a redder slope 
at longer wavelengths (see Fig. [6J. The differences illustrate the 
importance of K and M stars for the ISL. 

For the desired average scattered light spectrum, we run our 
scattering code for different combinations of zenithal optical 
depth, zenith distance, azimuth, and sidereal time for Rayleigh 
and Mie scattering (see Sect. 12.2b . The step sizes for the latter 
three parameters were 10°, 45°, and 2 h, respectively. Zenith dis- 
tances were only calculated up to 50°, since this results in a mean 
airmass of about 1 .25 for the covered solid angle, which is typi- 
cal of astronomical observations. The mean scattering intensities 
for each optical depth tq were then computed. This was trans- 
lated into a spectrum by using the relation between to and wave- 
length as provided by the Cerro Paranal extinction curve (see 
Sect. 12. II ). Finally, the mean ISL spectrum (normalised to the B- 
band flux) was multiplied. The sky model code does not change 
the final spectrum, except for the molecular absorption, which 
is adjusted depending on the bimonthly period (see Sect. 12. U . 
For the effective absorption airmass, the mean value of 1 .25 is 
assumed. 

The resulting spectrum shows an intensity of about 
13photons s _1 m~ 2 finT 1 arcsec~ 2 at 0.4jt/m. At 0.6yt/m, there is 
only half of this intensity (see Fig. [6]). The results are in good 
agreement with the findings of Bernstein et al. (120021 ) . 

3.3. Zodiacal light 

Zodiacal light is caused by scattered sunlight from interplanetary 
dust grains in the plane of the ecliptic. A strong contribution is 
found for low absolute values of ecliptic latitude B and heliocen- 
tric ecliptic longitude A - Aq. The brightness distribution pro- 
vided by Levasseur-Regourd & Dumont (1980) and Leinert et 
al. Q998 ) for 0.5 /mi shows a relatively smooth decrease for in- 
creasing elongation, i.e. angular separation of object and Sun. A 
striking exception is the local maximum of the so-called gegen- 
schein at the antisolar point in the ecliptic. The spectrum of the 
optical zodical light is similar to the solar spectrum (Colina et 
al. 119961 ), but slightly reddened (see Fig. [6}. We apply the rela- 
tions given in Leinert et al. (119981 ) to account for the reddening. 
The correction is larger for smaller elongations. Thermal emis- 
sion of interplanetary dust grains in the IR is neglected in the sky 
model, since the airglow components of atmospheric origin (see 
Fig. [TJ completely outshine it. In the optical, zodiacal light is a 
significant component of the sky model. A contribution of about 
50% is typical of the B and V bands when the Moon is down 
(see Sect. 15.21 ). At longer wavelengths, the fraction decreases 
due to the increasing importance of the airglow continuum (see 
Sect US. 

The model of the zodiacal light presented in Leinert et al. 
( 1998) describes the characteristics of this emission component 
outside of the Earth's atmosphere. Ground-based observations 
of the zodical light also have to take atmospheric extinction into 
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account. Since zodiacal light is an extended radiation source, the 
observed intensity is a combination of the extinguished top-of- 
atmosphere emission in the viewing direction and the intensity 
of light scattered into the line of sight (see Fig . |6j ■ The latter can 
be treated by scattering calculations as discussed in Sect. 12.21 
Scattering out of and into the line of sight leads to an effective 
extinction and this can be expressed by an effective optical depth 



T e ff - TQ.eff X - /ext T X 



(18) 



(cf. Eq.[TJ see also Bernstein et al. 2002). Consequently, the scat- 
tering properties of the zodiacal light can be described by the fac- 
tor / ext alone. This parameter shows only a weak dependence on 
To and, hence, wavelength. For Rayleigh scattering in the range 
from U to I band, we find an uncertainty of only about 4%. For 
Mie scattering from V to J, a similar variation is found. We ob- 
tain this result by calculating f ext for a grid of optical depths, 
zenith distances, azimuths, sideral times, and solar ecliptic lon- 
gitudes. We take the same grid as for scattered starlight (see 
Sect. 13.21 ) plus solar ecliptic longitudes in steps of 90°, i.e. for 
each of the four seasons, one data set was calculated. Scattering 
calculations were only carried out for solar zenith distances of 
at least 108°. This restriction excludes daytime and twilight con- 
ditions. As for scattered starlight, we consider target zenith dis- 
tances z up to 50° (see Sect. 13.21 ). For rare high zodiacal light 
intensities, we also let 50° < z < 70°. The z limits do not signifi- 
cantly affect /ext- A test with 0° < z < 70° showed a variation of 
/ext in the order of a few per cent only. 

For an efficient implementation in the ETC, we searched for 
a parametrisation of / ext that simplifies the treatment of scatter- 
ing and minimises the number of required input parameters. We 
found that the intial intensity of the zodiacal light in the viewing 
direction Iq indicates the tightest relation with / ext . The next-best 
parameter is the ecliptic latitude. Figure [7]shows / ext versus the 
logarithm of Iq as provided by Leinert et al. (19981 for Rayleigh 
and Mie scattering. Since the variation of / ext with to can be ne- 
glected, the optical depth is fixed to 0.27 for Rayleigh and 0.01 
for Mie scattering. These values correspond to the wavelengths 

0. 4 and 1 .2 jjm respectively, which represent typical wavelengths 
dominated by the two different scattering processes (see Fig. |2). 
The distribution of data points in Fig. Q exhibits a clear increase 
of /xt with /(). The slopes appear to be very similar for both scat- 
tering processes. The extinction reduction factor /xt,M for Mie 
scattering tends to be about 0.1 lower than / ex t.R for Rayleigh 
scattering. The difference increases with increasing Iq. For very 
low zodiacal light intensities in the viewing direction, / ext be- 
comes negative. Here, the scattering of light from directions of 
high zodiacal light brightness into the line of sight overcompen- 
sates the extinction. For viewing directions with higher Iq, this 
process is less effective and causes higher / ext . The investigated 
parameter space does not indicate values higher than about 0.75, 

1. e. the zodiacal light in the direction of astronomical observ- 
ing targets is always significantly less attenuated in the Earth's 
atmosphere than a point source (cf. Leinert et al. 19981. 

The good correlation of / ex t and Iq allows us to describe the 
reduction of the optical depth with Iq only. The unextinguished 
zodiacal light is derived from the ecliptic coordinates, which are 
given sky model parameters (see Table [TJ. Figure |7]displays the 
resulting linear fits to the data based on average / ext for 0. 1 dex 
bins of log /q. The parameters obtained by these split fits are 
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Fig. 7. Extinction reduction factor / ext for zodiacal light as a 
function of the line-of-sight top-of-atmosphere intensity Iq of 
the zodiacal light in dex. The intensities correpond to those 
given in Table 17 of Leinert et al. (119981 ). i.e. the unit is 
10~ 8 WirT^irT 1 sr~' (= 1.690photons s _1 rrr^rrT 1 arcsec~ 2 
at 0.5 fim). The factors were calculated for different ecliptic lon- 
gitudes of the Sun, sidereal times, and line-of-sight zenith dis- 
tances (up to 50° in general and up to 70° for log/o > 2.4) and 
azimuths. Only night-time data points with solar zenith distances 
greater than 108° are considered. For Rayleigh scattering (black 
'x' symbols), a zenithal optical depth of 0.27 is taken, corre- 
sponding to a wavelength of 0.4 yum for the Cerro Paranal ex- 
tinction curve (see Sect. 12. II ). For Mie scattering (green ' +' sym- 
bols), to is fixed to 0.01, a value reached at about 1.2/an. The 
figure also shows average / ext in 0. 1 dex Iq bins for both scatter- 
ing modes (Rayleigh: blue circles; Mie: red triangles). The re- 
sulting fits based on these average values are displayed by solid 
(Rayleigh) and dashed lines (Mie). 



for Rayleigh scattering and 



,/e> 



1.309 log/o 
0.468 log/o 



2.598 
0.702 



if log/o < 2.255, 
if log/o > 2.255 



(20) 



for Mie scattering with Iq in 10~ 8 Wm" 2 ^ 1 sr _1 . For the fits 
of the lower Iq, the / ext standard deviations of the data from the 
regression lines are 0.06 and 0.07 for Rayleigh and Mie scatter- 
ing respectively. For the higher /(>, systematic uncertainties are 
probably higher than statistical variations, since these values are 
rare when z < 70° and zq > 108°. 

The quality of the parametrisation in Eqs. [19] and [20] is il- 
lustrated by Fig. [8] which shows the average ratio of the fitted 
and originally calculated transmission for the entire parameter 
grid as derived from Eqs.[TJand[T8]in magnitudes. In almost the 
entire wavelength range, the accuracy of the mean is better than 
1%. Only at very short wavelengths, the systematic errors be- 
come higher. The intersection of the zero line at about 0.4/im is 
caused by the selection of to = 0.27 for the / ex t,R fitting. The lcr 
range suggests uncertainties of no more than 3 to 4% for most 
wavelengths. These results suggest that the parametrisation of 
T e ff does not cause higher errors than those expected from the 
uncertainties of the scattering calculations (see Sect. 12.21 ). The 
total errors from the simple treatment of multiple scattering and 
ground reflection, neglection of polarisation, simplified particle 
distributions, and uncertainties in the aerosol scattering param- 
eters should be in the order of 10% of / scat . This translates into 



/ext,R — 



1 .407 log Iq - 2.692 if log / < 2.244, 
0.5271og/ - 0.715 if log/ > 2.244 



(19) 
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Fig. 8. Ratio of the zodiacal light transmissions of the fits as 
shown in Fig. Q and the full scattering calculations in magni- 
tudes. The wavelength-dependent mean values and 1 <x devia- 
tions are displayed as solid and dashed lines respectively. 



uncertainties in the transmission of a few per cent at blue wave- 
lengths, which is somewhat higher than the uncertainties from 
the / ext fit. 

Finally, the reduction of the zodiacal light intensity by 
molecular absorption has to be considered. Strong absorption is 
only important for relatively long wavelengths, where scattering 
affects only a few per cent of the incoming light (see Fig. 
assuming zenith distances of typical astronomical observations. 
Hence, the direct component in the viewing direction is usu- 
ally much brighter than the scattered light from other directions. 
For this reason, we can safely assume that the target airmass is 
well suited to derive the transmission reduction by molecular 
absorption. Differences between target and effective scattering 
airmasses can be neglected. 



4. Airglow emission 

Atmospheric emission at wavelengths from the near-UV to the 
near-IR originates in the upper atmosphere, especially in the 
mesopause region at an altitude of about 90 km and in the iono- 
spheric F2 layer at about 270 km. In contrast to the thermal radia- 
tion in the mid-IR by atmospheric greenhouse gases in the lower 
atmosphere, this so-called airglow (see Khomich et al. 2008 for 
a comprehensive discussion) is caused by a highly non- 
emission process called chemiluminescence, i.e. by chemical re- 
actions that lead to light emission by the decay of excited elec- 
tronic states of reaction products. Consequently, atmospheric ra- 
diative transfer codes like LBLRTM (Clough et al. 120051 see also 
Sect. 12. Il l cannot be used to calculate this emission. Moreover, 
theoretical analysis is very challenging because of the high vari- 
ability of airglow on various time scales. The complex reactions 
are affected by the varying densities of the involved molecules, 
atoms, and ions in different excitation states. Therefore, airglow 
can best be treated by a semi-empirical approach. The section 
starts with a general discussion of airglow extinction (Sect. l4~Tl ). 
In Sect. 14.21 we present the data set that was used to derive the 
model. Detailed discussions of the line and continuum compo- 



Airglow 
layer 




Top of 
atmosphere 

Fig. 9. Illustration of the projected emission layer width used for 
the derivation of the beam-specific, incident airglow intensities 
in the 3D scattering calculations. The extension of the scatter- 
ing atmosphere is limited to the lower boundary of the airglow 
emission layer. Point O can be in a different plane from the other 
points (cf. Fig. 2]l. 



nents of the model can be found in Sects. 14.31 and 14741 respec- 
tively. 

4.1. Extinction of airglow emission 

The intensities of airglow lines and continuum rise with growing 
zenith distance by the increase of the projected emission layer 
width. This behaviour is expressed by the van Rhijn function 



K0) 



1 - 



R sin(z) 
R + h 



2\ 



(21) 



local thermodynamic equilibrium 



(van Rhijn |T92Tb . Here, z, R, and h are the zenith distance, the 
Earth's radius, and the height of the emitting layer above the 
Earth's surface (see Table |2), respectively. The airglow intensity 
is also affected by scattering and absorption in the lower atmo- 
sphere. Since the airglow emission is distributed over the entire 
sky, we can derive the effective extinction by means of scatter- 
ing calculations as described in Sect. 12.21 Some modifications 
have to be applied, since airglow emission originates in the at- 
mosphere itself and because of the van Rhijn effect. The top of 
atmosphere is reduced from the default H max - 200 km to 90 km, 
which is assumed to be representative of airglow emission lay- 
ers. The few atomic lines originating in the ionospheric F2 layer 
are neglected here (see Sect. 14.3) . The cut at the lower height is 
not critical, since the mass of the excluded range is negligible 
compared to the total mass of the atmosphere (see Eq. |9)- The 
projected emission layer width (see Fig. [9} is considered by cal- 
culating the zenith distance z ag of each incoming beam as seen 
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Fig. 10. Extinction reduction factor / ext for airglow as a func- 
tion of airmass (see Eg. 123} in dex. The factors were calculated 
for different line-of-sight zenith distances z (in 10° steps) and 
zenithal optical depths to. For Rayleigh (to < 0.99) and Mie 
scattering (to < 0.05), the values for to = 0.27 and 0.01 are 
marked by black 'x' and green '+' symbols respectively. For 
z < 60°, these values were used to derive airmass-dependent 
relations for / ext . The resulting fits are displayed by blue solid 
(Rayleigh) and red dashed lines (Mie). 
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Fig. 11. Ratio of the airglow transmissions of the fits as shown 
in Fig.QJJand the full scattering calculations in magnitudes. The 
wavelength-dependent results for line-of-sight zenith distances 
from 0° (black solid line) to 70° (red dashed line) in 10° steps 
are exhibited. 



from the corresponding entry point N at H msx . For our set-up, the 
resulting beam-specific airglow intensity can be approximated 
by 

/(Zag) = — j-t, (22) 

COS(z a g) 

since the centres of the scattering path elements S are always lo- 
cated well below the airglow emission layer height, which avoids 
z ag = 90°. The intensity Iq for a perpendicular incidence at N is 
assumed to be constant for the entire global airglow layer. The 



scattering calculations only depend on the zenithal optical depth 
and the target zenith distance. 

Like for the zodiacal light (see Sect. 13 . 3b , we can derive an 
optical depth reduction factor / ext (Eq. [TBI to describe the air- 
glow extinction properties. The representative optical depths are 
again 0.27 for Rayleigh scattering (6% uncertainty from U to 
/ band) and 0.01 for Mie scattering (4% variation from V to J 
band). The only remaining free parameter for / ext is the zenith 
distance or airmass. We take 



X ag = (l -0.972 sin 2 (z)) 



-0.5 



(23) 



(cf. Eq.fTTll. which corresponds to the van Rhijn formula (Eq.l2TI) 
for the assumed airglow layer height of 90 km. X ag is better 
suited for fitting / ext and Eq. [18] than the Rozenberg (1966) air- 
mass X (Eq. [2} used for the zodiacal light (see Sect. [33}, where 
the z dependence of / ext was negligible. FigureQJJshows the ex- 
tinction reduction factor as a function of the logarithm of X ig . 
There is a strong increase in / ext with airmass, and the factors 
for the Mie scattering tend to be about 0.15 lower than those for 
Rayleigh scattering. These trends are in good agreement with 
our findings for the zodiacal light. However, the typical val- 
ues of /ext are significantly lower. For most zenith distances 
relevant for astronomical observations, the factors are close to 
zero. Consequently, the effective extinction of airglow is very 
low (see Chamberlain 1961). For observations close to zenith, 
there is even a significant brightening of the airglow emission 
by the scattering into the line of sight from light with a long 
path length through the emission layer. Therefore, using a point- 
source transmission curve for airglow extinction correction will 
provide worse results than having no correction. As indicated by 
Fig. [10] a rough correction of airglow emission using the van 
Rhijn equation is reasonable up to zenith distances of 50° to 60°, 
depending on the ratio of Rayleigh and Mie scattering. At higher 
zenith distances, the compensation of extinction of direct light 
by scattered light from other directions becomes less efficient. 
At the largest z, the intensity can be more reduced by extinction 
than enhanced by the van Rhijn effect. For this reason, there is 
a zenith distance with a maximum airglow intensity, which de- 
pends on the optical depth and, hence, the wavelength. For the 
Cerro Paranal extinction curve (see Sect. 12.1) . we found a z max 
of about 70° at 0.35 yum and 85° at 0.7 pm. 

For the sky model, / ext was fit up to a zenith distance of 60°, 
where / ext shows an almost linear relation with the logarithm of 
X- dg (see Fig. [TO}. A worse agreement at large zenith distances is 
tolerable, since astronomical observations are rarely carried out 
at those angles. The parameters obtained by the fits are 



/ext,R = 1.669 log Xag- 0.146 
for Rayleigh scattering and 

/ ext , M = 1.732 logXa g - 0.318 



(24) 



(25) 



for Mie scattering. For the fixed optical depths and the fitted z 
range, the fit uncertainties of f ext are only about 0.01. Figure QT] 
displays the accuracy of the fits as a function of wavelength. For 
this plot, the Cerro Paranal extinction curve (see Sect. 12. U was 
used to convert optical depths into wavelengths. For the fitted 
range of zenith distances, the deviations between fitted and the 
true transmission curves are below 1%, at least for wavelengths 
longer than 0.4 /im. As expected, the deviations for z = 70° 
are significantly larger than for the smaller zenith distances. 
Nevertheless, for most prominent airglow lines (see Sect. 14.3b . 
the errors are below 2%, which also makes the parametrisation 
sufficient for large zenith distances. 
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Apart from scattering, molecular absorption in the lower at- 
mosphere affects the observed airglow intensity. However, in 
the optical these absorptions are usually small (see Fig. 
Exceptions are the airglow O 2 (b-X)(0-0) and O 2 (b-X)(l-0) 
bands at 762 and 688 nm, which suffer from heavy self- 
absorption at low altitudes, and can be observed as absorption 
bands only. As discussed in Sect. 13.31 the effective molecular 
absorption can be well approximated by taking the transmission 
spectrum for the target zenith distance. For the airglow contin- 
uum component, this is easily applied. For the line component, 
this is challenging, since airglow emission as well as telluric ab- 
sorption lines are very narrow. Therefore, an effective absorption 
for each airglow line has to be derived at very high resolution. 
Resolving the airglow lines requires resolutions in the order of 
10 6 . The absorption lines tend to be broader than the emission 
lines. LBLRTM (see Clough et al. 120051) easily produces spec- 
tra with the required resolution. For the airglow emission, we 
use a line list with good wavelength accuracy, as discussed in 
Sect. 14.31 In the upper atmosphere, broadening of lines by colli- 
sions can be neglected due to the low density (see Khomich et al. 
2008 ). Thus, the shape and width of airglow lines are determined 
by Doppler broadening, which can be derived by the molecular 
weight of the species and ambient temperature. The latter is in 
the order of 200 K in the mesopause region, which explains the 
low Doppler widths. For example, the FWHM of an OH line 
(see Sect. 14.31 ) at 0.8 pm is about 2picometers (pm). The few 
lines originating in the thermospheric ionosphere experience dis- 
tinctly higher temperatures of about 1000 K. The line absorption 
correction is applied as the multiplication of the initial line in- 
tensities by the suitable line transmissions from a pre-calculated 
list, which is derived from the library transmission spectrum cor- 
responding to the selected observing conditions (see Sect. 12.11 ). 

The airglow line absorptions show significant deviations 
from the continuum absorption at similar wavelengths. As al- 
ready mentioned, strong self-absorption is found for the O2 
ground state transitions in the optical. However, most other air- 
glow lines tend to be less absorbed than the continuum. In partic- 
ular, the H2O absorption bands (see Sect. 12. Q have little effect. 
This can be explained by the relatively low width of the telluric 
lines compared to the typical distance between strong absorption 
lines. Hence, the chance to find an airglow line at the centre of a 
strong absorption feature is relatively low. 

4.2. Data set for airglow analysis 
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The airglow analysis (see Sects. 1431 and !4.4l i and the sky model 
verification (see Sect. |5j were carried out using a sample of 
1186 FORS 1 long-slit sky spectra (see Fig. O- The FORS 1 
data were collected from the ESO archive and reduced by Patat 
(20080 The spectra were taken with the low/intermediate res- 
olution grisms 600B (12%), 600R (17%), and 300V, the latter 
with (57%) and without (14%) the order separation filter GG435. 
The set-ups cover the wavelength range from 0.365 to 0.89 fim. 
Wavelengths below 0.44 /mi are the least covered, since only 
600B and 300V spectra without GG435 can be used there. The 
FORS 1 spectra were taken between April 1999 and February 
2005, i.e. observations during a phase of low solar activity are 
not part of the present data set (see Fig. fT2l . The mean and 
standard deviation of the Penticton-Ottawa solar radio flux at 
10.7 crrQ (Covington 1 19691) are 153 and 35 solar flux units (sfu 



5 We excluded three spectra from the original sample of 1 189 spectra 
because of unreliable continua concerning strength and shape. 

6 ftp : // ftp . geolab . nrcan . gc . ca/data/ solar_f lux/ 
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Fig. 12. Date, time, and solar activity for the VLT FORS 1 obser- 
vations of the Patat (2008) spectroscopic data set. The solar ac- 
tivity is given by the solar radio flux measured at 10.7 cm in sfu. 
The FORS 1 spectra were taken with different instrument set- 
ups, which are indicated by different symbols and colours. The 
plot legend identifies the different set-ups by their covered wave- 
length ranges in fim. Open symbols indicate that the data were 
taken with the Moon above the horizon, whereas filled symbols 
refer to dark-time observations. 



= 0.01 MJy) respectively. The data set is also characterised by a 
mean zenith distance of 31°. The corresponding standard devi- 
ation is only 13° and the maximum zenith distance is 67°. The 
fraction of exposures with the Moon above the horizon amounts 
to 26%. 29 spectra were taken at sky positions with Moon dis- 
tances below 30° (see Sect. 13. II ). However, exposures affected by 
strong zodical light contribution are almost completely absent. 

The FORS 1 night-sky spectra were not extinction corrected 
and the response curves for the flux calibration were derived by 
the Cerro Tololo standard extinction curve (Stone & Baldwin 
119831 Baldwin & Stone 119841 We corrected all the spectra to 
the more recent and adequate Cerro Paranal extinction curve of 
Patat et al. d201 II see Sect. 12. U , assuming a mean airmass of 
1.25 for the spectrophotometric standard stars. This procedure 
resulted in corrections of +6% at the blue end to -1% at the red 
end of the covered wavelength range. 



4.3. Airglow lines 

In the following, we treat airglow lines in detail. At first, the 
origin of the airglow line emission, the classification of lines, 
and a derivation of a line list for the sky model are discussed 
(Sect. I43TTY In the second part, we investigate the variability of 
airglow lines (Sect. l43T2l . 



11 



S. Noll et al.: An atmospheric radiation model for Cerro Paranal 

Table 2. Basic properties of the Paranal airglow line and continuum (0.543 /mi) variability as derived from the Patat (20081 spectral 
data set. 



Pronprtv 


Unit 


TO 11 5577 


NalD 


roil 6300 6364 


OH(7 0W6 2) 


OiCh-XVO-H 


543 urn 

\J • ~J~^> LULL 


h layer 


km 


97 


92 
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94 


90 


N 

1 * spec 




1186 


1046 


1046 
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839 


874 


</>" 


phot s m arcsec 


3.6 


0.81 


3.5 


60. 


6.0 


80.* 






(1.6) 


(0.48) 


(3.4) 


(19.) 


(2.5) 


(31.*) 




R 


190. 


43. 


190. 


3200. 


320. 


4.3 rf 






(90.) 


(26.) 


(180.) 


(1000.) 


(130.) 


(1.6*) 


^sun 


sfu" 1 


0.0087 


0.0011 


0.0068 


0.0001 


0.0063 


0.0061 


p f 

1 season.min 




1 


1 


1 


5 


5 


1 


,/season,min 




0.86 


0.54 


0.30 


0.82 


0.71 


0.81 


p f 

1 season.max 




3 


3 


3 


6 


3 


3 


./season.max 




1.44 


1.84 


1.51 


1.14 


1.33 


1.24 



Notes. 

{a) Mean intensity and variability (in parentheses) for solar cycles 19 to 23. 

{b) In photons s~ m" 2 fim _l arcsec -2 . 

(c) 1 R (Rayleigh) « 0.018704 phot s" 1 nr 2 arcsec" 2 . 

W) inRnnr 1 . 

w Slope for solar activity correction relative to mean solar radio flux for cycles 19 to 23 (« 129 sfu). 

( -° Two-month bin with minimum or maximum correction factor: 1 = Dec/Jan, 2 = Feb/Mar, 3 = Apr/May, 4 = Jun/Jul, 5 = Aug/Sep, 6 = Oct/Nov. 
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Fig. 13. Variability classes for airglow emission lines. The fol- 
lowing groups are defined: (1) green OI, (2) NalD, (3) red OI, 
(4) OH, and (5) O2. The weak lines (green curves) are scaled by 
a factor of 30 for Na I D, red O I, and O2, and a factor of 10 for 
OH. 



4.3.1 . Line classes and line list 

The most prominent optical airglow line is [O I] 5577 (see 
Fig- El panel 1). Its dominant component at an altitude of about 
97 km is probably produced by the process described in Barth & 
Hildebrand (119611 1. The crucial reactions are 



O + O + M 



■M, 



OJ+O 



2 + 0( 1 S), 



and 



O('S) 0( 1 D) + 557.7 nm, 



(26) 



(27) 



(28) 



where M is an arbitrary reaction partner. The first reaction 
(Eq. [26} also plays an important role for the emission bands of 
molecular oxygen. In order to start the three body collision, oxy- 
gen molecules have to be split by hard UV photons: 



2 + hv -> O + O. 



(29) 



This process mainly occurs during the day by solar UV radia- 
tion, which implies that the oxygen airglow intensity depends 
on the solar activity and observing time. The O('D) state can 
lead to [O I] 6300 emission (see Fig. [13] panel 3). However, in 
the mesopause the deactivation of this metastable state is mostly 
caused by collisions. Hence, strong [O I] 6300 emission is re- 
stricted to the thermosphere at altitudes of about 270 km. In con- 
trast to the mesosphere, the excitation there is caused by disso- 
ciative recombination (Bates 1982), i.e. it depends on the elec- 
tron density: 

OJ + e - -> O + O('D). (30) 

The ionised molecular oxygen is produced by the charge transfer 
reaction 

+ + O^OJ+O. (31) 

The amount of ionised oxygen and free electrons depends on 
the solar radiation. The most striking features of the airglow are 
the Meinel OH bands (Meinel et al. |T954), which dominate at 
red optical wavelengths and beyond (see Fig. Q~3] panel 4). The 
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fundamental and overtone rotational-vibrational transitions are 
mainly produced by the Bates-Nicolet ( 1950) mechanism 

O3 + H -> OH* + 2 . (32) 

Most emission originates in a relatively thin layer at 87 km. 
Finally, the upper state of the prominent D lines of neutral 
sodium at 589.0 and 589.6 nm (see Fig. [13] panel 2), which orig- 
inate in a thin layer at 92 km, is probably excited by 

NaO + O -> Na( 2 P) + 2 (33) 

(Chapman |1939l ). NaO appears to be provided by 

Na + O3 -> NaO + 2 . (34) 

For studying airglow intensity and variability, we assigned 
the optical airglow lines to five classes, for which a similiar vari- 
ability is reasonable or could be proved by a correlation analy- 
sis (see Patat 2008). These classes are (1) green 01, (2) NalD, 
(3) red OI, (4) OH, and (5) 2 (see Fig. [13]). The first three 
classes correspond to atomic lines, the other two groups com- 
prise molecular bands. The O I lines are divided into two classes. 
The green OI line at 557.7 nm mainly originates at altitudes 
slightly below 100km, while the other significant OI lines (all 
being in the red wavelength range) tend to originate at altitudes 
greater than 200 km (see Khomich et al. 2008). Moreover, the 
kind of chemical reactions are also different, since the high al- 
titude emissions are usually related to reactions involving ions 
(see above). Another small group consists of the sodium D 
lines and a weak KI line at 769.9 nm (see Fig. Qj] panel 2). 
By far most of the airglow lines belong to the two classes of 
molecules which produce band structures by ro-vibrational tran- 
sitions. The OH electronic ground state (X-X) bands with up- 
per vibrational levels v' < 9 cover the wavelength range long- 
wards of about 0.5 //m with increasing band strength towards 
longer wavelengths. There are several electronic transitions with 
band systems for 2 in the optical. Relatively weak bands of 
the Herzberg I (A-X) and Chamberlain (A' -a) systems originate 
at near-UV and blue wavelengths (see Cosby et al. 2006]). At 
wavelengths longwards of 0.6 /mi, the atmospheric (b-X) band 
system produces several strong bands. However, those bands 
related to the vibrational v = level at the electronic ground 
state X are strongly absorbed in the lower atmosphere (see 
Sect. 14. 1) . Therefore, the only remaining strong band in the in- 
vestigated wavelength range is O 2 (b-X)(0-1) at about 864.5 nm 
(see Fig. [13] panel 5). Variability studies of 2 have to rely on 
the results from this band. 

For the identification of airglow lines and as an input line list 
for the sky model, we use the list of Cosby et al. (2006). It con- 
sists of 2805 entries with information on the line wavelengths, 
widths, fluxes, and line identities. It is based on the sky emis- 
sion line atlas of Hanuschik (2003 ), which was obtained from a 
total of 44 high-resolution (R ~ 45000) VLT UVES spectra. By 
combining different instrumental set-ups, the wavelength range 
from 0.314 to 1.043 /mi could be covered. The accuracy of the 
wavelength calibration is better than 1 pm (cf. Sect. 14.1) . The 
UVES spectra were flux calibrated by means of the Tug (11977) 
extinction curve. Since the use of a point-source extinction curve 
for the airglow leads to systematic errors (see Sect. 14.1) and the 
La Silla extinction curve was used for Cerro Paranal, we cor- 
rected the line extinction by means of the Patat et al. (201 ex- 
tinction curve and the recipes given in Eqs. l24l and l25l At the 
lower wavelength limit, this caused an extreme correction factor 
of 0.37 and in the range of the O 2 (b-X)(0-0) band a maximum 
factor of 4.5 due to the line molecular absorption correction (see 



Sect. 14.11 ). However, for the strongest 100 lines up to a wave- 
length of 0.92 /mi, we only obtain a mean correction factor of 
0.98. In a similar way as discussed in Sect. 14.21 for the Patat 
( 2008 ) data, we considered how a different extinction curve af- 
fects the response curves for flux calibration. Here, we obtain for 
the 100 brightest lines another correction of 0.97, which is rep- 
resentative of the factors for the entire wavelength range, which 
range from 0.95 to 1.00. Due to a gap at about 0.86 /mi, the 
Cosby et al. ([2006) line list misses a part of the important 2 (b- 
X)(0-1) band. There is an unpublished UVES 800U spectrum 
related to the study of Hanuschik (20031) that covers the gap. We 
used this spectrum to measure the line intensities. Determined by 
a few overlapping lines, the resulting intensities were then scaled 
to those in the line list. For the line identification, we considered 
2 and OH line data from the HITRAN database (see Rothman 
et al. 120091 ). For the final line list for the sky model, the inten- 
sities of each variability class were scaled to match the mean 
value derived from the variability analysis (see below). Finally, 
the line wavelengths were converted from air to vacuum by 

A vac = n /Ljr. (35) 

The refractive index 

0/ 2406030 15997 \ 

n - 1 + 1 ° ( 8342 - 13+ B0^ + 38^)' (36) 

where <x = A~ l and A is in /mi (Edlen [T966l ). This formula is 
also used internally in the sky model code if an output in air 
wavelengths is required (see Table[T}. 

4.3.2. Line variability 

In general, airglow lines show strong variability on time scales 
ranging from minutes to years. This behaviour can be explained 
by the solar activity cycle, seasonal changes in the temperature, 
pressure, and chemical composition of the emission layers, the 
day-night contrast, dynamical effects such as internal gravity 
waves and geomagnetic disturbances (see Khomich et al. 2008 ). 
In order to consider airglow variability in our sky model, we de- 
rived a semi-empirical model based on 1186 VLT FORS 1 sky 
spectra (Patat 2008; see Sect. 14.2) . The demands on an airglow 
variability model are twofold. First of all, it should include all 
major predictable variability properties. This excludes stochas- 
tic wave phenomena like gravity waves. Second, the derived 
parametrisation should be robust. For studying many variability 
triggers, about 1000 spectra is not statistically a high number. 
For this reason, only a few variables can be analysed. Since this 
neglection leads to higher uncertainties, an ideal set of param- 
eters has to be found that provides statistically significant, pre- 
dictable variations. For our airglow model, we studied the effect 
of solar activity, period of the year, and time of the night. The 
analysis was carried out for each of the five variability classes 
(see Fig. [13]). 

The first step for measuring the line and band fluxes was 
the subtraction of scattered moonlight, scattered starlight, and 
zodiacal light by using the recipes given in Sect. [3] The result- 
ing spectra were then corrected for airglow continuum extinc- 
tion as described in Sect. 14.11 For the flux measurements, con- 
tinuum windows were defined. Their central wavelengths are 
listed in Table [3] The widths were 4nm, except for 0.767 /mi, 
where 0.3 nm was chosen to reduce the contamination by strong 
OH lines and the O 2 (b-X)(0-0) band (see also Sects. @7T] and 
14.41 ). Then, the continuum fluxes were interpolated to obtain line 
and band intensities. The hydroxyl bands between 0.642 and 
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Fig. 14. Variability correction for the five airglow line classes and the airglow continuum (analysed at 0.543 pm) of the sky model. 
Upper panels: The variability is shown as a function of the bimonthly period (1 = Dec/Jan, 6 = Oct/Nov), time bin (third 
of the night), and solar activity measured by the solar radio flux. Lower panels: For bimonthly period and time bin, the relative 
uncertainties of the variability correction are displayed. 



0.858yum, i.e. OH(6-l) to OH(6-2), and the O 2 (b-X)(0-1) band 
between 0.858 and 0.872/vm were measured automatically. In 
a second iteration, this procedure was improved by subtracting 
the obtained O2 bands from the OH bands and vice versa be- 
fore the intensity measurement. This way, the contamination by 
undesired lines can be minimised. Apart from the bands, the 
lines [O I] 5577, NalD, [O I] 6300, and [O I] 6364 were anal- 
ysed. The continuum fit is very critical for these lines, because 
of the narrow width of the features and the significant contam- 
ination by OH lines. For this reason, the line measurements 
were also manually performed for a subsample of 15 spectra for 
each instrumental set-up. The results were then used to derive 
intensity-dependent correction functions for the automatic pro- 
cedure. Finally, the resulting band and line fluxes were corrected 
for the discrepancy between the molecular absorption of contin- 
uum and lines (see Sect. 14. U . These corrections were usually 
very small, since the strongest lines of each class do not overlap 
with significant telluric absorption. 

A variability model can be efficiently determined by deriv- 
ing a reference intensity and then applying multiplicative cor- 
rection factors for each variability parameter (cf. Khomich et al. 
2008). As a reference value for each variability class, we define 
the mean intensity at zenith for the five solar cycles 19 to 23, i.e. 
the years 1954 to 2007. The zenith intensity can be estimated by 
a correction for the van Rhijn effect (Eq.l2Tli. The correction fac- 
tors depend on the target zenith distance and the emission layer 



height (see Tabled. For describing the solar activity, we take the 
Penticton-Ottawa solar radio flux at 10.7 cm S 10.7 cm (Covington 
119691 see also Sect. 14.2b . For the defined period, we derive a 
mean flux of 129 sfu. We can also derive a mean S 10.7™ for the 
spectroscopic sample. For this purpose, we obtained the monthly 
S io.7cm averages corresponding to each spectrum. Due to the de- 
layed response of the Earth's atmosphere regarding solar activity 
of up to several weeks (see e.g. Patat 2008), monthly values are 
more suitable than diurnal ones. In the end, a mean S 10.7 cm of 
150 sfu could be derived for the entire FORS 1 data set. By per- 
forming a regression analysis for the relation between line inten- 
sity and S 10.7 cm for each line class, we were able to obtain cor- 
rection factors for a scaling of the lines to intensities typical of 
the reference solar radio flux. The resulting mean fluxes for the 
different classes were used to correct the fluxes in the reference 
line list (see Sect. 14.3. 11 1. To derive the scaling factors, we inte- 
grated over all the listed lines in a class that were measured in the 
spectra. This required adding up the OH intensities of all bands 
from (6,1) to (6,2) and summing up the fluxes of the [OI] 6300 
and [OI] 6364 lines, which have a fixed 3:1 intensity ratio. The 
extinction-corrected Cosby et al. (2006) intensities had to be cor- 
rected by factors between 0.43 for the red OI lines and 1.06 for 
NalD lines to reach the intensities for the reference conditions. 
These reference intensities are given in Table|2] They are slightly 
lower than those in Patat (2008) due to the lower mean solar ra- 
dio flux. 
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Fig. 15. Variability correction factors and their errors depending on bimonthly period (1 = Dec/Jan, 6 = Oct/Nov) and night time 
(third of the night: evening, midnight, morning) for the five airglow line classes and the airglow continuum (analysed at 0.543 yum) 
of the sky model. Data points only differing in the bimonthly period are connected by lines. 



The 5 10.7 cm-dependent linear fits to the airglow intensities 
(see also Patat 2008 ) are also used to parametrise a correction 
of solar activity for our sky model. Even though the parametri- 
sation is reasonable (cf. Khomich et al. 2008), there may be 
systematic deviations for monthly S 10.7 C m beyond the data set 
limits of 95 and 228 sfu (see Sect. 14.21 ). With the dependence 
on solar radio flux removed, seasonal and nocturnal variations 
can be investigated. Since the intensity has a complex depen- 
dence on these parameters (see Khomich et al. 2008), we di- 
vided the data set into 6 seasonal and 3 night-time bins. In the 
same way as described in Sect. 12. 1 I for the transmission curves, 
the two-month periods start with the December/January combi- 
nation. The night, as limited by the astronomical twilight, is al- 
ways divided into three periods of equal length. This means that 
in winter the periods are longer than in summer. For the latitude 
of Cerro Paranal (see Sect. [1}, the Dec/Jan night bins are about 
30% shorter than the Jun/Jul bins. Due to the additional con- 
sideration of specific averages over the entire year and/or night, 
the total number of bins defined for the sky model is 28 (see 
also TableQ]). Since the number of spectra in the individual bins 
differed by an order of magnitude, the average values over sev- 
eral bins could be biased. Therefore, all bins were filled with 
the same number of spectra by randomly adding spectra from 
the same bin, i.e. the data were cloned. This approach was also 
used for deriving the reference intensities. As a consequence, 
the mean solar radio flux of the data set changed from 153 (see 
Sect. 14.21 ) to 150 sfu. Finally, mean values and standard devia- 



tions were calculated for the individual and multiple bins using 
the original and bias-corrected data sets respectively. While the 
mean value provides an intensity correction factor for seasonal 
and/or night-time dependence, the standard deviation indicates 
the contribution of neglected variability components to the mea- 
sured intensity. It also depends on the data quality, uncertainties 
in the analysis, and the accuracy of other sky model components. 
Therefore, the standard deviation of each bin should be an upper 
limit for the unmodelled airglow variability in the investigated 
data set. 

The results of the airglow variability study are summarised 
in Table|2l and Figs. [T4l and [TBI The solar activity strongly influ- 
ences the intensities of OI and O2, where the slopes are 0.0063 
to 0.0087 sfu 1 . On the other hand, weak or no significant corre- 
lations are present for sodium and hydroxyl. These findings are 
in qualitative agreement with the relations given in Khomich et 
al. d2008l ). who provide 0.0025 to 0.0060 sfu 1 for the first and 
about 0.0015 sftr 1 for the second group from measurements at 
different observing sites (mainly in the Northern Hemisphere). 
The relatively low value of 0.0025 sfu 1 for the green OI line 
compared to 0.0087 sftr 1 from our model is not critical. A more 
precise latitude-dependent analysis of satellite-based data by Liu 
& Shepherd (120081 ) resulted in 0.0065 sftr 1 for the latitude range 
from -20° to -30°. For -30° to -40° they even obtain the same 
slope as in our model. The discrepancies in the results from lit- 
erature show that quantitative comparisons are difficult with dif- 
ferent observing sites, observing periods, instruments, and anal- 
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ysis methods. Solar activity directly affects the Earth's atmo- 
sphere by the amount of hard UV photons. The higher the pho- 
ton energy, the more the production rate depends on solar ac- 
tivity (see Nicolet 1989). This relation appears to be important 
for O I and O2 airglow emission, which depends on the dissocia- 
tion or ionisation of molecular oxygen by hard UV photons (see 
Sect.gTB. 

Significant dependence on the bimonthly period is found for 
all line classes (see Fig. [T4l>. NalD and the red OI lines indi- 
cate the strongest variability with a maximum intensity ratio of 
3.4 and 5.0 respectively, whereas OH shows the weakest depen- 
dence with a ratio of 1.4 (see Table[2]). Except for OH, the max- 
imum is reached in the autumn months April and May. For the 
atomic lines, the minimum arises in Dec/Jan, and for the molec- 
ular bands, it appears to be in Aug/Sep. The variations found re- 
flect the semi-annual oscilliations of airglow intensity (cf. Patat 
2008 ), which are characterised by two maxima and minima of 
different strength per year. Complex empirical relations based 
on long-term observations at many observing sites as given by 
Khomich et al. (2008) can be used to estimate this intensity vari- 
ation. However, since these relations depend on rough interpola- 
tions and extrapolations, they can only provide a coarse variabil- 
ity pattern for an arbitrary location on Earth. Nevertheless, these 
relations confirm the amplitude of the variations and the rough 
positions of the maxima and minima. For example, the peak-to- 
peak ratios for green OI, NalD, and OH are 1.4, 3.8, and 1.3, 
which agrees well with the sky model values 1.7, 3.4, and 1.4 
(see Table|2]). Apart from the obvious change of the incident so- 
lar radiation, the seasonal variations are related to changes in 
the layer heights of the atmospheric constituents, changes in the 
temperature profiles, and the dynamics of the upper atmosphere 
in general. 

Averaged over the year, the dynamical range of the night- 
time variations is in the order of 20 to 30%, which is distinctly 
smaller than the seasonal variations. Only the red O I lines show 
a clear decrease by a factor of 1.4 from the beginning to the 
end of the night. Usually, the diurnal variations are the largest in 
airglow intensity with possible amplitudes of several orders of 
magnitude (see Khomich et al. 120081) . Since we exclude daytime 
and twilight conditions with direct solar radiation, the variabil- 
ity is much smaller. In addition, the use of only three bins can 
smooth out the intensity variations. However, an investigation 
with the unbinned data (see also Patat 2008) and results from 
phot ometri c studies (Mattila et al. 119961 Benn & Ellison [19981 
Patat l2003l ) confirm the absence of strong trends. Averaging over 
all bimonthly periods appears to cancel part of the variability. 
This is illustrated in Fig. Q3] which shows the bimonthly de- 
pendence of the airglow intensity for each third of the night. 
At midnight, i.e. the maximum time distance to the twilight, 
the intensities of the mesopause lines tend to reach a minimum. 
However, the deviations from the intensities of the evening and 
morning bins appear to be a function of season. While in winter 
the differences reach a maximum, summer has the smallest dif- 
ferences. This behaviour can be explained by the night lengths 
at Cerro Paranal, which can vary by more than 30%. In sum- 
mer, the average length of time from the closest twilight is sig- 
nificantly smaller than in winter. This dependence supports the 
common treatment of seasonal and night-time variations in the 
sky model. 

The residual statistical variations of the airglow model, as 
depicted in Fig. [14] are mainly between 20 and 50% for the lines 
originating in the mesopause region. The smallest uncertainties 
are found for OH in summer, whereas the ionospheric O I lines 
indicate very strong unpredictable variations with an amplitude 
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Notes. The reference continuum intensity at 0.543 [im is 
79.8photons s _1 m 2 /jitT 1 arcsec~ 2 or4.27Rmrr'. 

comparable to the intensity itself (see also Patat 2008). This be- 
haviour is caused by the aurora-like intensity fluctuations from 
geomagnetic disturbances. This effect is particularly important 
for Cerro Paranal, because of its location close to one of the two 
active regions about 20° on either side of the geomagnetic equa- 
tor (Roach & Gordon [T973l 

4.4. Airglow continuum 

The airglow continuum is the least understood emission compo- 
nent of the night sky. In the optical, the best-documented process 
is a chemiluminescent reaction of nitric oxide and atomic oxy- 
gen in the mesopause region as proposed by Krassovsky (119511 ): 

NO + O -> N0 2 + hv (37) 

N0 2 + O -> NO + 2 . (38) 

This reaction appears to dominate the visual and is expected to 
have a broad maximum at about 0.6 yum (see Sternberg & Ingham 
[19721 von Savigny et al. |1999l Khomich et al. [2008b . Reactions 
of NO and ozone could be important for the airglow continuum 
at red to near-IR wavelengths (Clough & Thrush [T967I Kenner 
& Ogryzlo 1984). In particular, the reaction 

NO + 0| -> NOJ + 2 (39) 

is expected to significantly contribute to the optical by contin- 
uum emission with a broad maximum at about 0.85 /urn (Kenner 
& Ogryzlo [HHIl Khomich et al. l2008l . Finally, a pseudo con- 
tinuum by molecular band emission produced by 

Fe + 3 -> FeO* + 2 (40) 

(West & Broida [T975l ) could be an important component of the 
continuum between 0.55 and 0.65 /um (Jenniskens et al. 2000; 
Evans et al. 120101 Saran et al. l20TT1 ). 

The airglow continuum was analysed in a similar way as the 
airglow emission lines. As reported in Sect. 14.3.21 the deriva- 
tion of line intensities already required the measurement of air- 
glow continuum fluxes in suitable windows (see Table[3]). These 
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Fig. 16. Airglow/residual continuum (solid line) and its variabil- 
ity (dashed line) relative to the reference wavelength 0.543 /im. 



results were refined by subtracting possible contributions from 
airglow lines by applying the line emission model as described 
in Sect. 14.31 For the variability analysis, the Patat (2008) data 
set (see Sect. 14.21 1 was restricted to the 874 FORS 1 spectra, 
where the Moon was below the horizon. This avoids large er- 
rors in the continuum flux from the subtraction of a bright and 
uncertain component (see Sect. 13. 11 1. Even though the absolute 
uncertainties of the other components are significantly smaller, 
the derived continuum is a residual continuum. It is affected 
by uncertainties in zodical light brightness, scattered starlight, 
airglow line emission, molecular absorption, possible uncon- 
sidered minor components, scattered light from dispersive ele- 
mental and flux calibration of the spectra. While the other sky 
model components should cause an uncertainty in the order of a 
few per cent only (for 0.767 /mi this could be somewhat higher, 
see Sect. 14.3.2b . the uncertainty in the flux calibration could be 
significantly more. Since theoretical components are subtracted 
from the observed spectra to obtain the airglow continuum, er- 
rors in the flux calibration have a larger effect on the resulting 
continuum than on the total spectra. For example, for a 50% con- 
tribution of the continuum to the total flux, a flux error of 10% 
would convert into a continuum error of 20%. 

For the subsequent analysis, we defined a reference con- 
tinuum wavelength. Due to its availability in all spectroscopic 
modes and the absence of contaminating emission lines, we se- 
lected 0.543 /im. In the following, we studied the continuum rel- 
ative to this wavelength, i.e. the mean continuum was derived 
by scaling the fluxes to the reference solar radio flux of 129 sfu 
(see Sect. 14.3.2b by using only the linear regression results and 
binning corrections derived for 0.543 ji/m. This approach means 
that a fixed shape of the continuum is assumed. This simplified 
variability model can be more easily used in the sky model. As 
shown by the last column of Table[3] there are relatively low vari- 
ations of the continuum flux relative to the reference flux for the 
investigated data set. For most wavelengths, the standard devia- 



7 This effect can cause line profiles with wide Lorentzian wings (see 
Ellis & Bland-Hawthorn 20081. However, we did not detect significant 
wings in the profiles of the strong airglow lines or lines in the wave- 
length calibration lamp spectra. In view of the relatively low contrast 
between lines and continuum for our spectra (cf. the near-IR regime), 
we can conclude that a possible continuum contamination of more than 
a few per cent is unlikely, even at the reddest wavelengths. 



tion divided by the mean flux 07 / (/) is lower than 20% and only 
reaches 30% at the reddest wavelengths. In view of all the uncer- 
tainties discussed above, which also contribute to these numbers, 
the assumption of a fixed continuum slope is sufficient for an ap- 
plication for the ETC sky background model. This assumption is 
also supported by the high correlation factors for different con- 
tinuum windows measured by Patat (2008). It should be noted 
that the number of available spectra decreases towards the mar- 
gins in the analysed spectral range due to the different instru- 
mental set-ups (see Sect. 14.21 1, In particular, wavelengths short- 
wards of 0.44 yum have only 239 spectra, i.e. 27% of the sample 
used for 0.543 yum. The change in the sample size as well as dif- 
ferent systematic errors in the response curves for the different 
set-ups cause systematic uncertainties in the airglow continuum 
shape. By comparing the mean fluxes in overlapping wavelength 
regions of the different set-ups, we estimate an uncertainty in the 
order of 10%. 

We obtained a mean airglow continuum for Cerro Paranal 
with an intensity of about 4.3 R nirT 1 at 0.543 yum (for R[ayleigh] 
units see Table 0. The corresponding variability is 1.6Rnm _1 . 
Considering the latitude of Cerro Paranal, this result is in good 
agreement with observations that suggest an increase of the 
airglow continuum from 1 - 2RnnT 1 at the equator to 10 - 
20 R nm" 1 in the polar region for similar wavelengths (see Davis 
& Smith [T9651 Khomich et al. [2008b . As illustrated in Fig. [16] 
and listed in Table [3] the mean airglow continuum is charac- 
terised by a broad bump at about 0.6 fim and an increase in in- 
tensity towards longer wavelengths with a steep slope at about 
0.8 //m. The continuum in the bump region appears to be mainly 
produced by two processes: a relatively smooth NO + O con- 
tinuum (see Eq. l37l > and a more structured Fe + O3 continuum 
(see Eq. 1401 showing a peak at about 0.6/im (see Evans et al. 
1201 Oi l. If the bump, which indicates a width of about 0.1 //m, 
could be completely explained by the iron oxide reaction, this 
component would contribute up to 1 /3 of the mean airglow con- 
tinuum flux (see also Sect. 15. U . Towards longer wavelengths, it 
is not clear whether the observed increase in intensity is mainly 
caused by relatively narrow peaks or a smooth increase (see 
Sternberg & Ingham [T972l Sobolev [T978l Content [19961 . Up to 
0.872yum, our measurements suggest a monotonic increase be- 
yond the 0.6/i/m bump. Intensities of about 15 Rnrrr 1 , as found 
by Sternberg & Ingham (|1972| l and Sobolev (1978) at wave- 
lengths longwards of 0.8 yum, are in good agreement with our 
mean spectrum, with the same intensity at about 0.86//m. Also, 
the intensities of Sternberg & Ingham (11972) and of the sky 
model mean airglow continuum are very similar at our refer- 
ence wavelength in the V band. A possible production meach- 
anism for the strong increase at about 0.8 fim could be a reac- 
tion of nitric oxide with excited ozone (see Eq. [39l . for which 
a peak wavelength of about 0.85 yum is expected (see Kenner & 
Ogryzlo 1984; Khomich et al. 120081) . Then, the strong intensity 
increase would be the blue wing of the bump from this reaction. 
The slope does tend to decrease where the peak is predicted, and 
the continuum variations there relative to 0.543 yum (see 07 / (/) 
in Table [3) are the highest. This could be another argument for 
different chemical processes producing the emissions at 0.55 and 
0.85 yum. 

As in the case of the airglow lines, the intensity of the sky 
model airglow continuum measured at 0.543 yum depends on 
zenith distance, solar radio flux, period of the year, and period 
of the night. The results are shown in Table [2] and Figs. [T4l and 
[TBI The dependence of the continuum emission on solar radio 
flux is significant and comparable to O I (green and red) and O2. 
Since OH does not indicate a significant dependence on S 10.7 cm. 
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Fig. 17. Comparison of the sky model (upper panel: thick black 
line) and an observed FORS 1 300V spectrum (thin red line) with 
moderate lunar contribution (see TableQ}. The uncertainty of the 
sky model due to airglow variability is also shown (green). The 
lower panel exhibits the ratio of the sky model and the observed 
spectrum. 



the continuum must be the main source for the dependence of 
the broad-band fluxes on solar activity (cf. Walker 1988 Patat 
2003 2008). Concerning the seasonal and nocturnal variations, 
the continuum at 0.543 ^m tends to show the best agreement 
with OH, i.e. the variability of these parameters is relatively low. 
The maximum ratio of the bimonthly mean intensities is only 
about 1.5 (cf. PatatHH. 



5. Discussion 

After describing all the optical model components in the previ- 
ous sections, we now discuss the quality of the sky model as a 
whole by comparing it to the Patat ( 2008) data set of sky spectra 
and results found in the literature. First, sky model spectra and 
observed spectra are compared (Sect. 15.1) . Second, a system- 
atic quality evaluation by means of broad-band photometry of 
our model with data, and a comparison to the night-sky bright- 
nesses from the literature are discussed in Sect. 15.21 Finally, we 
evaluate the improvements achieved by our sky model compared 
to the current ESO ETC sky background model in the optical 
(Sect.O- 

5.1. Comparison of spectra 

In order to evaluate the quality of our sky model, spectra were 
computed for the different instrumental set-ups and observing 
conditions of the Patat (2008) data set (see Sect.l4~2l). A detailed 
description of the sky model input parameters can be found in 
Table [TJ Before the full data set is discussed, we show a few 
typical sky model spectra in comparison with the corresponding 
observed spectra (Figs.[T7lto[T9l. 
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Fig. 18. Comparison of the sky model and an observed FORS 1 
300V spectrum with strong lunar contribution. For an explana- 
tion of the different curves, see Fig. [171 
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Fig. 19. Comparison of the sky model and an observed FORS 1 
300V spectrum without lunar contribution, but strong airglow 
emission lines. For an explanation of the different curves, see 

Hg.na 



In the optical, the relevant sky model components are scat- 
tered moonlight, scattered starlight, zodiacal light, and airglow 
lines and continuum of the upper atmosphere (see Fig. [TJ. 
Figure[T7]compares a typical FORS 1 spectrum exhibiting mod- 
erate Moon contribution to the corresponding sky model. This 
and the FORS 1 spectra of the other figures were taken with the 
300V grism without order separation filter. This set-up provides 
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Fig. 20. Deviations between the sky model and the observed 
FORS 1 spectra in magnitudes. The mean magnitude differ- 
ence (upper panel), the standard deviation (middle panel), and 
the wavelength-dependent number of considered spectra (lower 
panel) are shown for 1 nmbins. Results for the full spectroscopic 
data set (red) and for spectra with the Moon below the horizon 
(black) are displayed. 



the maximum wavelength range. Due to the missing order sepa- 
ration filter, there can be a slight (5= 10%) contamination by the 
second order spectrum at red wavelengths. However, this does 
not affect the quality of the sky model (see Sect. 15.2b . because 
of the small number of spectra used in the analysis. Figure [T7] 
demonstrates that the different sky model components agree rea- 
sonably well with the observed spectrum. Apart from the red- 
dest part of the spectrum, deviations of 10 to 20% are typical. 
The best agreement is found for the continuum in the visual. 
Residuals of strong sky lines in the ratio plot in the lower panel 
appear more prominent than the real values due to slight differ- 
ences in the instrumental profiles of the model and the observed 
spectrum. 

Figure [18] shows the effect of a strong Moon contri- 
bution on the sky model. Scattered moonlight dominates 
the continuum at all wavelengths in the displayed range. 
At blue wavelengths, the example spectrum reaches up to 



2000photons s m /vm arcsec . In contrast, during moon- 
less nights the values are about one order of magnitude lower 
(see Fig.[T9b. Even though the deviations between the sky model 
and the observed spectrum are not larger than 10% for most of 
the continuum, there is a systematic trend. For a wide wave- 



length range, the sky model appears to be too bright by a rel- 
atively constant factor. Apart from deviations of the Moon illu- 
minance from the model in Eq. [15] this discrepancy could be 
explained by uncertainties in the Mie scattering model, since 
this changes slowly with wavelength (see Fig. [2). Moreover, the 
angular distance between Moon and target p = 44° is small 
enough for significant Mie scattering. Similar trends with pos- 
itive as well as negative deviations are observed in other cases 
of high lunar sky brightness. As already mentioned in Sect. 13.11 
the amplitude of the differences tends to increase for decreasing 
p. Nevertheless, the uncertainties of the extended Krisciunas & 
Schaefer (1991) model are usually not critical for spectra with 
low to moderate lunar contribution, as it is illustrated in Fig.fTTl 

A spectrum without scattered moonlight is shown in Fig. |T9| 
In this case, the contiuum is mainly a combination of zodia- 
cal light and airglow continuum. Scattered starlight only adds 
a few per cent to the total continuum emission. In the B and 
V bands, the continuum of dark-time exposures is usually well 
reproduced by the sky model. At longer wavelengths, the uncer- 
tainties increase from the significant (and partly unpredictable) 
variability of the airglow continuum and the simplifying assump- 
tion of a fixed continuum shape (see Sect. 14.41 ). The three spectra 
shown indicate deviations between 20 to 40% at about 0.82 pm. 
By chance, the sky model is brighter than the observed spectrum 
in all cases. Airglow emission lines can vary significantly, de- 
pending on the line class, observing direction, solar activity, sea- 
son, and time of night. Even though there is always a stochastic 
component, which cannot be modelled, the sky model accounts 
for such variations quite well. The uncertainty of the line mod- 
elling for the spectra shown is in the order of 20%. The spec- 
tra in Figs. [TTI and [T9l were taken at the same solar radio flux 
(205 sfu) and in the last third of the night, i.e. the conspicu- 
ous differences between both spectra can mainly be explained 
by seasonal variations. This is confirmed by the fact that the sky 
model can roughly follow the striking changes in the intensities 
of the red O I lines (cf. FM and JJ periods for red O I morning 
intensities in Fig. [15]). 

For a more systematic comparison between the model and 
the data, we computed the mean and standard deviations for the 
entire sample and a subset with the Moon above the horizon. The 
number of spectra involved depends on the wavelength and does 
not exceed 1 186 and 874 respectively. The results are displayed 
in Fig.[20]in magnitudes. For both samples, the systematic model 
deviations (upper panel) are close to zero. This is very satisfying, 
but not unexpected. The airglow line and continuum components 
were derived based on the FORS 1 spectra, which are also used 
in the comparison. The lack of clear spectral features in the upper 
panel of Fig. [20] suggests that the model quality is similar at all 
wavelengths. As already discussed in Sect. 14.41 the real accuracy 
of the model also depends on the quality of the flux calibration 
of the Patat (2008) data. The uncertainties can be in the order 
of 10%, and so this appears to be the main source of systematic 
errors in the model intensities. 

The standard deviation of the magnitude differences between 
the model and observed data cr^ m , shown in the middle panel 
of Fig. [20] indicates an interesting wavelength dependence. For 
the dark-time sample, the continuum uncertainty ranges between 
0.16 and 0.24 mag. Including observations with scattered moon- 
light increases these values to 0. 19 and 0.28 mag. The additional 
variation is mainly caused by a few spectra with large uncertain- 
ties in the Mie scattering of moonlight (see above). These results 
suggest that the relative strong continuum deviations at 0.82 pm, 
as in Figs. [17] and [19] are relatively high compared to the entire 
sample. The spectrum of cr^ m shows a relatively broad feature 
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at about 0.6 /im, which is reminiscent of the shape of the FeO 
pseudo continuum (see Jennikens et al. 2000; Evans et al. 2010; 
Saran et al. 1201 11 1. This result supports the assumption that the 
airglow continuum bump at 0.6 ytzm (see Fig. [TBI is mainly pro- 
duced by iron oxide (see Sect. 14.4b and varies differently than 
the underlying NO2 continuum, for which the intensitiy changes 
at 0.543 /im are considered by the sky model. Fig. [20] indicates 
a FeO-related variability of about 15%. Since the mean bump 
contributes about 20% of the total flux at 0.6 //m, the variation 
is almost comparable to it. These strong variations are in agree- 
ment with previous results (see Evans et al. 120101 Saran et al. 
1201 11 1. A minor source of continuum variability at 0.6 fj.m is the 
Chappuis ozone bands (see Fig. |2|, which indicate a variation 
of 0.02 mag airmass -1 (see Patat et al. 2011). The middle panel 
of Fig. [20] does not exhibit strong signatures of OH bands (see 
Fig. [T3] panel 4). The bands are clearly detectable only in the 
range of low continuum uncertainty. Consequently, the quality 
of the OH modelling is comparable or even better than the con- 
tinuum modelling. This is particularly important for the qual- 
ity of the full sky model at long wavelengths, where OH emis- 
sion dominates the spectra. The sky model is less able to model 
atomic lines, as shown in Fig. [20] The maximum variability of 
0.7 mag is found at [O I] 6300. In addition, other lines of the ther- 
mospheric ionosphere, which are usually difficult to recognise 
(see Fig. [13] panel 3), significantly contribute to <x Am . For exam- 
ple, NI at 520 nm and OI at 616, 656, and 777 nm are clearly 
detectable. This confirms the conclusion of Sect. 14.3.21 stating 
that the intensities of ionospheric lines are difficult to predict 
due to their sensitivity to geomagnetic disturbances (see Roach 
& Gordon lT973l . 

5.2. Comparison of photometry 

To perform a more quantitative comparison of our sky model and 
observed data, and for an easier comparison to other studies, we 
derived magnitudes from broad-band filter fluxes. Specifically, 
we used the standard photometric system consisting of U, B, V, 
R, and / (see Bessell [19901 . For the Patat ([2008 ) data and their 
corresponding sky models, only a grism-dependent subset of fil- 
ters can be taken. The spectrum needs to cover 90% of the 
filter curve. An exception is the U filter, which extends below 
0.365 yum, the lower wavelength limit of the bluest spectral set- 
up (see Sect. 14.21 1. Therefore, the results for this filter are ques- 
tionable. 

Figure [2JJ and Table |4] show the results of the computations 
for the full sample of FORS 1 spectra and a subsample with the 
Moon below the horizon (cf. Sect. I5.lt . Only a small fraction 
(~ 1/4) of the spectra contribute to the mean values and stan- 
dard deviations of the U and B filters. Most spectra cover the 
range of the V, R, and / filters. For both samples and all filters, 
the mean magnitude differences are small (^ 0.05 mag). Even 
though the sky model tends to be slightly fainter at the blue end 
and slightly brighter at the red end, there are no serious sys- 
tematic deviations from the observed spectra (cf. Sect. 15. lb . The 
standard deviations cr Am are relatively similar for the different 
filters. For the full sample, they range from 0.19 mag for the V- 
band to 0.24 mag for the / band. Consequently, the accuracy of 
the sky model is in the order of about 20%. These values can be 
compared to the magnitude variations in the measured data <x b s , 
which correspond to cr Am for an optimum sky model of time- 
invariant flux. For the full sample, the deviations are distinctly 
larger. They range from 0.37 to 0.56 mag. These large deviations 
are mainly due to the observations with a strong moonlight con- 
tribution (see Sect. 13. U . For dark-time conditions, (Am) ranges 
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Fig. 21. Histograms for the deviations of the sky model from the 
FORS 1 data in mag for the filters U, B, V, R, and /. The filled 
histograms show only data with the Moon above the horizon. 



from 0. 17 mag for B to 0.21 mag for U, and cr Am indicates differ- 
ences between 0. 18 mag for B and 0.29 mag for /. In particular, 
for the bands V to /, a time-invariant sky model would be signifi- 
cantly worse than our sky model (see also Sect. 15.3) . This shows 
the improvements from a variable zodiacal light, airglow lines, 
and airglow continuum model. 

Finally, we compare the sky model to literature data. 
Magnitudes in the five broad-band filters U to / were calculated 
using standardised observing conditions. Specifically, zenith, 
New Moon, ecliptic pole, annual average, and mean solar ac- 
tivity, i.e. 130 sfu (see Sect. 14.3.21 . were assumed. The resulting 
magnitudes are provided in Table[5] For a better comparison with 
published data, the table also contains sky model results for 90 
and 180 sfu. The reference data in Table [5]originate from Mattila 
et al. ([19961 for La Silla, Walker (fl987T > and K risciunas et al. 
([20071 for Cerro Tolol o, Benn & Ellison ([T9981 for La Palma, 
and Patat (2003, 2008) for Cerro Paranal. In general, there is a 
good agreement between the sky model and observed sky bright- 
nesses. By comparing magnitudes with similar average solar ra- 
dio fluxes, we find typical deviations in the order of 0.1 mag. 
An exception is the U filter with typical differences of about 
0.3 mag, which could be explained by the lower quality of the 
sky model at near-UV wavelengths due to the lack of data (see 
Sect. 14. 21 ). Also, it is not clear to what extent components like zo- 
diacal light or scattered starlight, which significantly affect short 
wavelengths, contribute to the observed data. The corresponding 
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Table 4. Comparison of the sky model and the Walker ( 1987 ) magnitudes with the Patat (120081) FORS 1 data. 



Filter All No Moon 





N 


(Am) 


C Am 


C obs 


<Am w ) 




N 


(Am) 


O" Am 


fobs 


(A/W w > 


f Am,w 


u 


303 


+0.053 


0.230 


0.521 


-0.113 


0.725 


239 


+0.042 


0.213 


0.214 


-0.059 


0.703 


B 


303 


+0.046 


0.205 


0.561 


+0.269 


0.430 


239 


+0.048 


0.174 


0.182 


+0.196 


0.407 


V 


979 


+0.000 


0.192 


0.563 


+0.356 


0.409 


727 


+0.006 


0.175 


0.255 


+0.244 


0.322 


R 


1046 


+0.009 


0.199 


0.447 


+0.169 


0.365 


762 


-0.009 


0.177 


0.260 


+0.083 


0.287 


I 


839 


-0.003 


0.238 


0.365 


+0.413 


0.326 


615 


-0.040 


0.208 


0.287 


+0.365 


0.294 



Notes. Mean values (Am) and standard deviations cr A „, for magnitude differences between the model and observed data (Am = m raodel - m obs ) 
and standard deviations of the magnitudes of the observed data cr b s are shown. Moreover, mean values (Am w ) and standard deviations cr Am w for 
magnitude differences between the Moon phase related data of Walker and the Patat data are displayed. Results are listed for the full and dark-time 
samples. The number of spectra involved is indicated by N. 

Table 5. Typical sky model and literature night-sky brightnesses in mag arcsec~ 2 for zenith, no Moon, faint zodiacal light, and 
different solar radio fluxes. 



Source 


Site 


S 10.7 cm" 


U 


B 


V 


R 


/ 


Benn& Ellison (1998J 


La Palma 


80 


22.0 


22.7 


21.9 


21.0 


20.0 


Walker ilWf) 


Cerro Tololo 


90 


22.0 


22.7 


21.8 


20.9 


19.9 


Krisciunas et al. ( 2007 1 


Cerro Tololo 


130 


22.1 


22.8 


21.8 


21.2 


19.9 


Mattila et al. (119961 


La Silla 


150 




22.8 


21.7 


20.8 


19.5 


Patat J2008T> 


Cerro Paranal 


160 


22.4 


22.7 


21.7 


20.9 


19.7 


Patat (20031 


Cerro Paranal 


180 


22.3 


22.6 


21.6 


20.9 


19.7 


Sky model* 


Cerro Paranal 


90 


22.3 


22.9 


22.0 


21.2 


19.8 






130 


22.1 


22.8 


21.8 


21.0 


19.7 






180 


21.9 


22.6 


21.6 


20.9 


19.6 



Notes. 

{a] Solar radio flux measured at 10.7 cm in sfu. 
<b) See TablefTJfor model parameters. 



intensities can differ from those of the sky model for the assumed 
ideal conditions. The published photometry appears to confirm 
the lack of reproducibility of [/-band magnitudes. Unexpectedly, 
the sky brightness seems to increase with decreasing solar activ- 
ity. It should be noted that photometric samples are often too 
small to cover all of the airglow variability, which causes addi- 
tional uncertainties. For the airglow continuum, which is a major 
contribution to the night-sky brightness in the optical, the inten- 
sity significantly depends on the geographic latitude (see Davis 
& Smith 1965 ). In summary, we can assume that our sky model 
calculates reliable night-sky brightnesses, at least for the levels 
of solar activity that could be analysed. 

A comparison of the sky model with other studies can also 
be enlightening in terms of the night-sky variability, particularly 
the striking dependence of the sky brightness on solar activ- 
ity (Walker [T9M1 Mattila et al. [E>96l Krisciunas [19971 Patat 
2003 120081 ). For the range of photometric fluxes in UVI un- 
der dark-time conditions at zenith, Patat (2008 ) found 0.9, 1.3, 
and 1.7 mag. For the same filters, he estimated a flux increase of 
0.6, 0.3, and 0.2 mag for a change in the solar radio flux from 80 
to 240 sfu. In contrast, Walker (119881 ) gives AB « 0.8 mag and 
AV =s 1 .0 mag for a similar solar activity change. About 0.5 mag 
are reported by Mattila et al. (1996 ) for blue wavelengths and 
by Krisciunas et al. d 19971) for the V band. For a typical bright- 
ness of the zodiacal light of 1.7 times the intensity at the ecliptic 
pole and the same range of solar radio fluxes as used by Patat 
([20081 . the sky model gives 0.7, 0.5, 0.6, 0.6, and 0.3 mag for the 
UBVRI filters. These values are consistent with the results from 
the other studies. The significant discrepancy of 0.3 mag to the 
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Patat ( 2008 ) result for the V filter might be explained by the rela- 
tive small sample of only 148 images and the general difficulties 
in separating the different sky radiation components by means of 
photometric data. For B to R, the solar activity appears to be re- 
sponsible for about half the night-sky variability. Since the zodi- 
acal light variations cause a 1 cr uncertainty of about 0. 1 mag, as 
estimated from the FORS 1 data set, seasonal, nocturnal, and un- 
predictable short-term airglow variations should generate most 
of the remaining variability in V and the adjacent bands. 

5.3. Comparison to current ETC sky model 

The main motivation for this study is to achieve a significant im- 
provement over the current ESO ETC for a more efficient use of 
the VLT time. As discussed in Sect. [T] the current sky emission 
model for optical instruments consists of the Moon phase de- 
pendent, photometric sky brightness table of Walker ( 1987 ) and 
an instrument-specific template spectrum scaled to the flux in 
the photometric filter closest to the centre of a given spectral set- 
ting. It does not consider any change in the relative contributions 
of the different sky emission components, which can be quite 
strong as discussed in the previous sections. Thus, we can ex- 
pect that the model continuum level significantly deviates from 
the true one, even if the corresponding filter flux is correct. Our 
new sky model with variable model components is an important 
improvement. 

For a more quantitative comparison of the current and our 
new sky model, we derived the Walker (119871 ) magnitudes for 
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the observing conditions of the Patat (2008 ) FORS 1 data set (see 
Sect. 14.2b . The resulting U to I magnitudes can be compared to 
the measured ones to analyse the quality of the two models (see 
also Sect. l5.21 >. Table[4]exhibits the mean deviation of the Walker 
model (Am w ) and the standard deviation of the difference o"A m , w 
for the full and the dark-time sample. Since the Walker magni- 
tudes depend on the Moon phase, one can expect that the model 
based on his data set would have a lower scatter than an opti- 
mum time-invariant model given by <x obs . Except for the U band, 
which is difficult to compare (see Sect. 15.2b . this is true. For the 
full data set, the standard deviations are lower by 0.04 (/) to 
0.15 mag (V). However, this improvement is small compared to 
the corresponding results for our sky model. For the B and V 
bands cr Am w is more than twice as large as <x Am . Moreover, there 
are significant systematic mean deviations of 0.1 to 0.4 mag. For 
our sky model, the maximum value is 0.05 mag. Even though 
this benefits from the use of the same data for the derivation of 
the airglow line and continuum model, the offset with the Walker 
model is striking. It can be explained by the fact that the Moon 
phase is the only variable parameter. Therefore, the systematic 
changes of other parameters can cause significant magnitude off- 
sets. For example, the fainter magnitudes of the Walker model 
can partly be explained by the very low solar activity (about 
90sfu), whereas the Patat (2008 ) data set has a mean solar ra- 
dio flux of 153 sfu. 

The dark-time data provided by Table |4] show a decrease of 
o"Am,w by 0.02 (B) to 0.08 mag (R) compared to the values for 
the full sample. Excluding the U band, the typical scatter of the 
Walker model is about 0.3 mag. This is worse than the results for 
our sky model and even a time-invariant model. For the B band, 
the difference is more than 0.2 mag. For the / band, the situation 
is better, since the Walker model almost reaches the standard 
deviation of the time-invariant model. Nevertheless, cr Am of our 
sky model is still about 0. 1 mag lower. The unsatisfactory per- 
formance of the Walker model is caused by the unconsidered, 
varying amount of scattered moonlight during a Moon phase. 
Even for phases close to Full Moon, dark-time observations are 
possible. For this reason, the scatter of the Walker model is larger 
than for a time-invariant model. The difference depends on the 
average contribution of scattered moonlight to the sky bright- 
ness, which is higher at shorter wavelengths. This effect causes 
an overestimation of the sky brightness for dark-time observa- 
tions, which is indicated by a comparison of (Am w ) for the full 
and dark-time samples. 

In summary, we can state that our sky model performs better 
than the current ESO ETC sky model for the optical by about 
several tenths of a magnitude. The quality difference should be 
even higher for a single wavelength instead of a broad-band fil- 
ter. Even if the input parameters of our new sky model (see 
Table [TJ could not optimally be constrained in the planning 
phase, due to uncertainties in the scheduling of the observations, 
there would be a significant improvement for the prediction of 
the signal-to-noise ratio for astronomical observations. This jus- 
tifies the effort in developing an advanced sky model. 



6. Summary and outlook 

In this paper, we presented a sky model for the ESO observing 
site at Cerro Paranal to help predict the effects of the Earth's 
atmosphere on astronomical observations in the optical. This 
sky model includes atmospheric extinction, scattered moonlight, 
scattered starlight, zodiacal light, and airglow emission lines and 
continuum. 



- The atmospheric extinction is characterised by three compo- 
nents: Rayleigh scattering, aerosol extinction, and molecular 
absorption. For the former, we use the formula given by Liou 
(2002). The aerosol extinction is based on the parametri- 
sation of Patat et al. d201 II) for Cerro Paranal. The molec- 
ular absorption is calculated by the radiative transfer code 
LBLRTM for atmospheric profiles optimised for the observ- 
ing site. 

- Scattered moonlight is derived from the Krisciunas & 
Schaefer (1991) model for the V band, which is extrapolated 
for the entire optical range using the solar spectrum and the 
Cerro Paranal extinction curve. 

- For scattered starlight, a mean spectrum was derived by 3D 
single scattering calculations with multiple scattering correc- 
tions, using the Toller ( 119811 ) integrated starlight distribution 
and the mean spectrum of Mattila ( 198DJ. 

- The zodiacal light is computed based on the model given by 
Leinert et al. ( 11998I ). Scattering of the zodiacal light in the 
Earth's atmosphere is also treated with 3D scattering calcu- 
lations. 

- Airglow emission lines and continuum are highly variable. 
We used a set of 1 186 VET FORS 1 sky spectra (Patat |2008) 
taken during a wide range of conditions to parametrise a 
model for five line and one continuum variability classes. 
The line wavelengths and intensities are based on the atlas 
of Hanuschik (2003 ), improved by the results of the variabil- 
ity analysis. For the airglow scattering, 3D calculations were 
used as well. 

After summing over all the components, we then compared 
the sky model with observed sky spectra. We found that our 
model is accurate to about 20%. This could only be achieved 
by having variable model components, especially the airglow 
model, which is a significant improvement over previous mod- 
els. The mean sky brightnesses and variations derived from our 
model are in good agreement with results from other photomet- 
ric studies. Compared to the current ETC sky model, the accu- 
racy is increased by several tenths of a magnitude. 

There are new results from this sky model study: 

- We found scattering scaling relations for the extended radia- 
tion sources zodiacal light and airglow. These are new sim- 
ple relations to derive effective extinction curves. The cor- 
responding correction factors for the point-source extinction 
curve only depend on the unextinguished line-of-sight inten- 
sity of the radiation source. 

- We characterised the variability of the airglow line and con- 
tinuum emission at Cerro Paranal based on several input pa- 
rameters, such as period of the year, time of night, and solar 
radio flux. 

- We derived long-term mean airglow line and continuum in- 
tensities. For the latter, we benefit from the modelling of 
the intensity of the other sky model components. We could 
determine the shape of the optical airglow continuum with 
high accuracy. For the controversial wavelength range long- 
wards of 0.7 /im, we find a smooth intensity increase. This 
could be produced by chemical reactions different from those 
dominating the shorter wavelengths. The well-known peak at 
0.6/mi was identified as strongly varying. 

The model presented in this paper is aimed at the ETC for 
optical VET imagers and spectrographs. In the near/mid-IR, ther- 
mal emission from the telescope, airglow lines (mainly OH and 
O2 bands) and continuum and, in particular, emission/absorption 
bands from greenhouse gases, originating mainly in the lower at- 
mosphere, are dominating the sky background of ground-based 
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observations. The scattering-related components become negli- 
gible. The change in significance between the components as 
well as the necessity to use IR data sets for the analysis requires 
different modelling techniques to achieve reliable ETC predic- 
tions in the IR. The IR sky radiation model will be described in 
a forthcoming paper. 

The sky model can also be used as a basis for procedures 
aiming at reconstructing data by removing the sky signature. 
This purpose requires the adaption of the sky model to enable 
the fitting of observed data. Preliminary studies show promising 
results for this kind of application, which will also be presented 
in future papers. 
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